Larvae Biomineralization Block 1
Load libraries
Load data Name the data you export from FileMaker Pro by their exact table names and save them as CSVs, e.g. Larvae_Morphology.csv
#setwd("~/Repos/LarvaeTransGen2018")
setwd("~/R/GitHub/LarvaeTransGen2018/data") #Elise's working directory
#upload all of the data tables for the larvae experiment
Larvae_Morphology <- read.csv("../data/Larvae_Morphology.csv") #contains data from CellProfiler from the larvae outlines
Barcode_Jar <- read.csv("../data/Barcode_Jar.csv") #contains data on CrossID, seatable, and whether or not larvae were present
Block_ID <- read.csv("../data/Block_ID.csv") #contains data on the block
Cross<- read.csv("../data/Cross.csv") #contains data on the female and male IDs for the cross and whether the cross was for QG or Meth
Fert_QG<- read.csv("../data/Fert_QG.csv") #contains data on fertilization counts, the JarIDs for the crosses
Header_WC<- read.csv("../data/Header_WC.csv") #contains data on the header tanks used to fill the jars
Larvae_Calcification <- read.csv("../data/Larvae_Calcification.csv") #contains data on filter weights and counts
Larvae_Counts <- read.csv("../data/Larvae_Counts.csv") #contains data on sedgewick rafter larvae counts
Larvae_WC <- read.csv("../data/Larvae_WC.csv") #contains data on larvae jar water chemistry
Parent_ID <- read.csv("../data/Parent_ID.csv") #contains data on the adults at the time of shucking. Includes Parent ID assignments
Larvae_Cilia<- read.csv("../data/Larvae_Cilia.csv") #contains data on cilia extrusion of the larvae photographed for morphology
Egg_Morphology<- read.csv("../data/Egg_Morphology.csv") #contains data on egg sizes from CellProfiler
Adult_Sample<- read.csv("../data/Adult_Sample.csv") #contains data on adult oysters at the time of collection
WC_Standard<- read.csv("../data/WC_Standard.csv") #contains data on the standards used for water chemistry
Tank_WC<- read.csv("../data/Tank_WC.csv") #contains data on the adult tank water chemistryMake fields factors:
Adult_Sample$AccTankID<- as.factor(Adult_Sample$AccTankID)
Adult_Sample$TrtTankID<- as.factor(Adult_Sample$TrtTankID)
Barcode_Jar$JarSeatable<- as.factor(Barcode_Jar$JarSeatable)
Barcode_Jar<- subset(Barcode_Jar, JarSeatable =="2" | JarSeatable =="4" | JarSeatable =="6")
Tank_WC$ParentTrt<- as.factor(Tank_WC$ParentTrt)
Tank_WC$Tank<- as.factor(Tank_WC$Tank)
Tank_WC$Folder<- as.factor(Tank_WC$Folder)
WC_Standard$CRM<- as.factor(WC_Standard$CRM)
Larvae_WC$DateFilter<- as.factor(Larvae_WC$DateFilter)Calculate average larval jar water chemistry based on the subset of bottles that were run on the VINDTA.
#start by joining Barcode_Jar and Larvae_WC datasets
JarChem<- merge(x=Barcode_Jar, y= Larvae_WC, by="JarID", all.x=TRUE)
#combine Date filter and JarTrt columns to get a unique combo value for each
JarChem<-unite(JarChem,BlockTrt, c("DateFilter","JarTrt"), sep='', remove=F)
#get only jars with larvae
JarChem<- subset(JarChem, Larvae =="TRUE")
JarChem<- subset(JarChem, DateFilter =="20180713")
boxplot(JarpHSW~BlockTrt, data=JarChem, ylim=c(6,8))
#get means for the parameters
JarMeans<- ddply(JarChem, .(BlockTrt), numcolwise(mean, na.rm=T))
#remove the unnecessary columns from JarMeans
JarMeans<- subset(JarMeans, select= c(BlockTrt,JarAlkCalc, JarpCO2, JarHCO3, JarCO3, JarCO2, JarOmegaCalcite, JarOmegaArg))
#add Est for "estimated" to the column names
colnames(JarMeans)<- paste("Est", colnames(JarMeans), sep="")
#add columns to JarChem for mean carb values
JarChem<- merge(x=JarChem, y=JarMeans, by.x="BlockTrt", by.y="EstBlockTrt", all.x=TRUE)
#use ifelse statement to fill in the na values in the original carb fields
JarChem$JarpCO2<- ifelse((is.na(JarChem$JarpCO2)), paste(JarChem$EstJarpCO2), JarChem$JarpCO2)
JarChem$JarAlkCalc<- ifelse((is.na(JarChem$JarAlkCalc)), paste(JarChem$EstJarAlkCalc), JarChem$JarAlkCalc)
JarChem$JarHCO3<- ifelse((is.na(JarChem$JarHCO3)), paste(JarChem$EstJarHCO3), JarChem$JarHCO3)
JarChem$JarCO3<- ifelse((is.na(JarChem$JarCO3)), paste(JarChem$EstJarCO3), JarChem$JarCO3)
JarChem$JarCO2<- ifelse((is.na(JarChem$JarCO2)), paste(JarChem$EstJarCO2), JarChem$JarCO2)
JarChem$JarOmegaCalcite<- ifelse((is.na(JarChem$JarOmegaCalcite)), paste(JarChem$EstJarOmegaCalcite), JarChem$JarOmegaCalcite)
JarChem$JarOmegaArg<- ifelse((is.na(JarChem$JarOmegaArg)), paste(JarChem$EstJarOmegaArg), JarChem$JarOmegaArg)
#get the carb parameters to numeric
JarChem$JarpCO2<- as.numeric(JarChem$JarpCO2)
JarChem$JarAlkCalc<- as.numeric(JarChem$JarAlkCalc)
JarChem$JarHCO3<- as.numeric(JarChem$JarHCO3)
JarChem$JarCO3<- as.numeric(JarChem$JarCO3)
JarChem$JarCO2<- as.numeric(JarChem$JarCO2)
JarChem$JarOmegaCalcite<- as.numeric(JarChem$JarOmegaCalcite)
JarChem$JarOmegaArg<- as.numeric(JarChem$JarOmegaArg)Look at the water chemistry data for adult tanks
#subset the tank data to only include B1
Tank_WCB1<- subset(Tank_WC, Block == "B1")
Tank_WCB1$Tank<-droplevels(Tank_WCB1)$Tank #drop unused levels
#test to see if tanks are significantly different.
boxplot(WCa_pHDIC~Tank, data=Tank_WCB1, na.omit=TRUE)
#check for significant differences in saturation state between tanks within treatments
Ca1<- lmer(WCa_pHDIC~ParentTrt + (1|Tank), data=Tank_WCB1)## boundary (singular) fit: see ?isSingularCa2<- lmer(WCa_pHDIC~ 1 + (1|Tank), data= Tank_WCB1)
anova(Ca1, Ca2) #select Ca1## refitting model(s) with ML (instead of REML)## Data: Tank_WCB1
## Models:
## Ca2: WCa_pHDIC ~ 1 + (1 | Tank)
## Ca1: WCa_pHDIC ~ ParentTrt + (1 | Tank)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## Ca2 3 9.9301 14.327 -1.9651 3.9301
## Ca1 4 -21.1940 -15.331 14.5970 -29.1940 33.124 1 8.646e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Ca3<- lm(WCa_pHDIC~ParentTrt, data=Tank_WCB1)
anova(Ca1, Ca3) #select Ca3## refitting model(s) with ML (instead of REML)## Data: Tank_WCB1
## Models:
## Ca3: WCa_pHDIC ~ ParentTrt
## Ca1: WCa_pHDIC ~ ParentTrt + (1 | Tank)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## Ca3 3 -23.194 -18.797 14.597 -29.194
## Ca1 4 -21.194 -15.331 14.597 -29.194 0 1 1#check assumptions
par(mfcol=c(2,2))
plot(Ca3) #seems to meet assumptions of homoscedasticity and linearity. Normality could look better. 
par(mfcol=c(1,1))
acf(resid(Ca3)) #this looks good. 
summary(Ca3) #There is a significant effect of ParentTrt, but not tank according to best model. ##
## Call:
## lm(formula = WCa_pHDIC ~ ParentTrt, data = Tank_WCB1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.38809 -0.04330 -0.00918 0.06594 0.25937
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.20034 0.03959 30.32 < 2e-16 ***
## ParentTrt2800 -0.85889 0.05599 -15.34 9.66e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1584 on 30 degrees of freedom
## Multiple R-squared: 0.8869, Adjusted R-squared: 0.8832
## F-statistic: 235.3 on 1 and 30 DF, p-value: 9.659e-16#try a transformation to see if it helps normality
Ca1log<- lmer(log(WCa_pHDIC)~ParentTrt+(1|Tank), data=Tank_WCB1)
Ca3log<- lm(log(WCa_pHDIC)~ParentTrt, data=Tank_WCB1)
anova(Ca1log, Ca3log) #still choose Ca3log## refitting model(s) with ML (instead of REML)## Data: Tank_WCB1
## Models:
## Ca3log: log(WCa_pHDIC) ~ ParentTrt
## Ca1log: log(WCa_pHDIC) ~ ParentTrt + (1 | Tank)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## Ca3log 3 -15.952 -11.5545 10.976 -21.952
## Ca1log 4 -13.952 -8.0888 10.976 -21.952 0 1 1par(mfcol=c(2,2))
plot(Ca3log) #now meets assumptions 
par(mfcol=c(1,1))
acf(resid(Ca3log))#looks good
summary(Ca3log)##
## Call:
## lm(formula = log(WCa_pHDIC) ~ ParentTrt, data = Tank_WCB1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.37383 -0.08255 0.00641 0.14423 0.29896
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.16588 0.04434 3.741 0.000773 ***
## ParentTrt2800 -1.25239 0.06270 -19.974 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1773 on 30 degrees of freedom
## Multiple R-squared: 0.9301, Adjusted R-squared: 0.9277
## F-statistic: 399 on 1 and 30 DF, p-value: < 2.2e-16#create dataframe that has the mean values for each tank
TankMeans<- ddply(Tank_WCB1, .(Tank, Block, ParentTrt), numcolwise(mean, na.rm=T))Join the data into dataframe for Egg analysis
#First dataframe will be for examining egg morphology for all eggs traced; multiple values for each adult
ParentInfo<- merge(x=Adult_Sample, y=Parent_ID, by="AdultID", all.x=TRUE)
ParentInfo<- merge(x=ParentInfo, y=TankMeans, by.x="TrtTankID", by.y="Tank", all.x=TRUE)
Eggs<- merge(x=Egg_Morphology, y=ParentInfo, by.x="FemaleID", by.y="ParentBlockID", all.x=TRUE)
Eggs$FemaleID<- as.factor(Eggs$FemaleID)
#subset the data to only include B1
EggsB1<- subset(Eggs, Block =="B1")Visualize the data for eggs to start
#get the original range of the eggs
range(EggsB1$EggDiamum)## [1] 15.35819 92.72958#look at the egg measurements to see outliers
boxplot(EggDiamum~FemaleID, data=EggsB1)
boxplot(EggDiamum~ParentTrt, data=EggsB1)
Remove egg outliers
#check for outliers using code from https://www.r-bloggers.com/identify-describe-plot-and-remove-the-outliers-from-the-dataset/. The function uses Tukey's method to ID ouliers ranged above and below the 1.5*IQR.
outlierKD <- function(dt, var) {
var_name <- eval(substitute(var),eval(dt))
na1 <- sum(is.na(var_name))
m1 <- mean(var_name, na.rm = T)
par(mfrow=c(2, 2), oma=c(0,0,3,0))
boxplot(var_name, main="With outliers")
hist(var_name, main="With outliers", xlab=NA, ylab=NA)
outlier <- boxplot.stats(var_name)$out
mo <- mean(outlier)
var_name <- ifelse(var_name %in% outlier, NA, var_name)
boxplot(var_name, main="Without outliers")
hist(var_name, main="Without outliers", xlab=NA, ylab=NA)
title("Outlier Check", outer=TRUE)
na2 <- sum(is.na(var_name))
cat("Outliers identified:", na2 - na1, "n")
cat("Propotion (%) of outliers:", round((na2 - na1) / sum(!is.na(var_name))*100, 1), "n")
cat("Mean of the outliers:", round(mo, 2), "n")
m2 <- mean(var_name, na.rm = T)
cat("Mean without removing outliers:", round(m1, 2), "n")
cat("Mean if we remove outliers:", round(m2, 2), "n")
response <- readline(prompt="Do you want to remove outliers and to replace with NA? [yes/no]: ")
if(response == "y" | response == "yes"){
dt[as.character(substitute(var))] <- invisible(var_name)
assign(as.character(as.list(match.call())$dt), dt, envir = .GlobalEnv)
cat("Outliers successfully removed", "n")
return(invisible(dt))
} else{
dt[as.character(substitute(var))] <- invisible(var_name)
assign(as.character(as.list(match.call())$dt), dt, envir = .GlobalEnv)
cat("Outliers successfully removed", "n")
return(invisible(dt))
}
}
outlierKD(EggsB1, EggDiamum) #I opted to remove outliers for now. 
## Outliers identified: 48 nPropotion (%) of outliers: 11.5 nMean of the outliers: 57.63 nMean without removing outliers: 62.4 nMean if we remove outliers: 62.95 nDo you want to remove outliers and to replace with NA? [yes/no]:
## Outliers successfully removed n#Outliers identified: 48 nPropotion (%) of outliers: 11.5 nMean of the outliers: 57.63 nMean without removing outliers: 62.4 nMean if we remove outliers: 62.95 n
yes #to remove outliers from EggsB1## Error in eval(expr, envir, enclos): object 'yes' not found#count how many eggs are left per female after removing the outliers.
outrem<- aggregate(EggDiamum~FemaleID, data=EggsB1, FUN=length)
#EF07_B1 had 10 eggs that were included the rest were removed. Keep EF07_B1 in the analysisVisualize eggs without outliers
#Plot again without outliers
par(mfcol=c(1,1))
boxplot(EggDiamum~FemaleID, data=EggsB1) 
boxplot(EggDiamum~ParentTrt, data=EggsB1)
#get the range now
range(EggsB1$EggDiamum, na.rm=TRUE)## [1] 49.95329 75.96955#note that the eggs were filtered through a 70 um filter, but it appears that larger ones made it through.We are going to keep these in for now because they aren't that outrageously large.
#plot histograms of the egg data with outliers removed
ggplot(EggsB1, aes(x=EggDiamum, color=ParentTrt))+
geom_histogram(alpha=0.5, position="identity")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 48 rows containing non-finite values (stat_bin).
eggmeans<- aggregate(EggDiamum~ParentTrt, data=EggsB1, FUN=mean) #400 = 63.46233; 2800= 62.15902
eggmedians<- aggregate(EggDiamum~ParentTrt, data=EggsB1, FUN=median)#400 = 62.91200; 2800= 61.80109
view(eggmeans)
view(eggmedians)
#I will use the mean egg size. Make datasets without the egg size outliers
meanegg<- ddply(EggsB1, .(FemaleID), numcolwise(mean, na.rm=T)) #get mean egg size
#select only columns that we need for egg morph
meanegg<- subset(meanegg, select=c("FemaleID", "EggDiamum"))
FemAdults<- merge(x=ParentInfo, y= meanegg, by.x="ParentBlockID", by.y= "FemaleID", all.y=TRUE, all.x=TRUE)
#Next dataframe will be for one line per traced larva, includes all the data for that jar, including average egg size for that jar
Larvae<-merge(x=Larvae_Morphology, y=JarChem, by="JarID", all.x=TRUE)
#use gather function in tidyr to make the Fert_QG Jar ID columns in long form vs. wide form
Fert_QG_long <- Fert_QG %>% tidyr::gather(Key, JarID, JarID1:JarID6)
#use the gathered data frame to merge with Larvae
Larvae<-merge(x=Larvae, y=Fert_QG_long, by.x= c("JarID", "CrossID"), by.y=c("JarID", "CrossID"), all.x=TRUE) #did the join by two columns so CrossID column wasn't duplicated
Larvae<- merge(x=Larvae, y=Cross, by="CrossID", all.x=TRUE)
Larvae<- merge(x=Larvae, y=FemAdults, by.x= "FemaleID", by.y="ParentBlockID", all.x=TRUE)
Larvae<- merge(x=Larvae, y= Larvae_Counts, by="JarID", all.x=TRUE)
#Remove jars without larvae
Larvae<- subset(Larvae, Larvae=="TRUE")
#get the growth per day for the larva
Larvae$GrowthPerDay<- (Larvae$LarvaeDiamum-Larvae$EggDiamum)/3
#Final dataframe will be one line per Jar, includes all data for the jar, including the average larva size per jar and cilia
LarvByJar<-merge(x=Larvae_Calcification, y=Larvae_Cilia, by="JarID", all.x=TRUE, all.y=TRUE)
LarvByJar<- merge(x=LarvByJar, y=JarChem, by="JarID", all.x=TRUE)
LarvByJar<- merge(x=LarvByJar, y=Fert_QG_long, by.x=c("JarID", "CrossID"), by.y=c("JarID","CrossID"), all.x=TRUE)
LarvByJar<- merge(x=LarvByJar, y=Cross, by="CrossID", all.x=TRUE)
LarvByJar<- merge(x=LarvByJar, y=FemAdults, by.x="FemaleID", by.y="ParentBlockID", all.x=TRUE)
LarvByJar<- merge(x=LarvByJar, y=Larvae_Counts, by="JarID", all.x=TRUE)
#get the means for larvae morphology and add to LarvByJar
MeanLarvMorph<- ddply(Larvae_Morphology, .(JarID), numcolwise(mean, na.rm=T))
SELarvMorph<- ddply(Larvae_Morphology, .(JarID), numcolwise(se, na.rm=T)) #get standard error of larva size
MeanLarvMorph<- full_join(MeanLarvMorph,SELarvMorph, by=c("JarID", "JarID"),suffix=c("","SE")) #join the mean and se larvae data
LarvByJar<- merge(x=LarvByJar, y=MeanLarvMorph, by="JarID", all.x=TRUE)
#Now let's add the survival data. See notes on "Larvae Survival" for more info on how we decided to count total larvae in the jar
#Make columns for the two larvae counts for 1_5 and 2
LarvByJar$v1_5SurvCount<- LarvByJar$SWRTotal+LarvByJar$F2LarvaeCount
LarvByJar$v2SurvCount<- LarvByJar$F1LarvaeCount+LarvByJar$F2LarvaeCount
#use ifelse statement to define the LarvaeSurvived column
LarvByJar$LarvaeSurvived<- ifelse(LarvByJar$ProtocolVersion =="1", paste(LarvByJar$TotalLarvae), ifelse(LarvByJar$ProtocolVersion =="1_5", paste(LarvByJar$v1_5SurvCount), paste(LarvByJar$v2SurvCount)))
LarvByJar$LarvaeSurvived<- as.numeric(LarvByJar$LarvaeSurvived)## Warning: NAs introduced by coercion#get the total larvae for the control jars and exposed jars for each family and then get the ratio of exposed to control to get an idea of how good a male and female pair are at larvae surviving. Can also use this to create reaction norms.
SurvDat<- ddply(LarvByJar, .(CrossID, ParentTrt, JarTrt, MaleID, FemaleID, TrtTankID, AccTankID, Block), numcolwise(mean, na.rm=TRUE))
SESurvDat<- ddply(LarvByJar, .(CrossID, ParentTrt, JarTrt, Block), numcolwise(se, na.rm=TRUE)) #get the standard errors for each mean count
SESurvDat<- subset(SESurvDat, select= c(CrossID, JarTrt, LarvaeSurvived)) #rename the columns
colnames(SESurvDat)<- paste("SE", colnames(SESurvDat), sep="")
#merge the standard error data to the SurvDat dataframe
SurvDat<- merge(SurvDat, SESurvDat, by.x= c("CrossID", "JarTrt"), by.y= c("SECrossID", "SEJarTrt"))
#only keep the relevant fields because otherwise this is a huge dataframe
SurvDat<- subset(SurvDat, select= c(CrossID, JarTrt, ParentTrt, MaleID, FemaleID, TrtTankID, AccTankID, LarvaeSurvived, Block, SELarvaeSurvived, WCa_pHDIC))
SurvDat<- subset(SurvDat, Block=="B1")
SurvDat<- subset(SurvDat, CrossID!="")
#use the spread function to go from long form to wide form for the survival data
SurvDatToWide<- unite(SurvDat, Value, c(LarvaeSurvived, SELarvaeSurvived), sep="_")
SurvWideDat<- spread(SurvDatToWide, JarTrt, value= Value)
#now split the mean and SE apart again
SurvWideDat<- separate(SurvWideDat, Control, into= c("MeanLarvSurvivedCon", "SELarvSurvivedCon"), sep="_")
SurvWideDat<- separate(SurvWideDat, Exposed, into= c("MeanLarvSurvivedExp", "SELarvSurvivedExp"), sep="_")
#make the means and ses numeric
SurvWideDat$MeanLarvSurvivedExp<- as.numeric(SurvWideDat$MeanLarvSurvivedExp)
SurvWideDat$MeanLarvSurvivedCon<- as.numeric(SurvWideDat$MeanLarvSurvivedCon)
SurvWideDat$SELarvSurvivedExp<- as.numeric(SurvWideDat$SELarvSurvivedExp)## Warning: NAs introduced by coercionSurvWideDat$SELarvSurvivedCon<- as.numeric(SurvWideDat$SELarvSurvivedCon)## Warning: NAs introduced by coercionSurvWideDat$RatSurv<- SurvWideDat$MeanLarvSurvivedExp/SurvWideDat$MeanLarvSurvivedCon
SurvWideDat$SurvChange<- (SurvWideDat$MeanLarvSurvivedCon-SurvWideDat$MeanLarvSurvivedExp)/SurvWideDat$MeanLarvSurvivedCon
SurvWideDat$CrossID<- as.factor(SurvWideDat$CrossID)
SurvWideDat$FemaleID<- as.factor(SurvWideDat$FemaleID)
SurvWideDat$MaleID<- as.factor(SurvWideDat$MaleID)
LarvByJar<- subset(LarvByJar, Larvae=="TRUE")
#now remove the data we do not want to analyze.
LarvaeDat<- subset(Larvae, Block == "B1")
LarvaeDat$JarID<- as.factor(LarvaeDat$JarID)
LarvaeDat$AccTankID<- as.factor(LarvaeDat$AccTankID)
LarvaeDat$TrtTankID<- as.factor(LarvaeDat$TrtTankID)
LarvaeDat$JarSeatable<- as.factor(LarvaeDat$JarSeatable)
LarvaeDat$JarTrt<- as.factor(LarvaeDat$JarTrt)
LarvaeDat$ParentTrt<- as.factor(LarvaeDat$ParentTrt)
LarvaeDat$FemaleID<- as.factor(LarvaeDat$FemaleID)
LarvaeDat$MaleID<- as.factor(LarvaeDat$MaleID)
LarvaeDat$CrossID<- as.factor(LarvaeDat$CrossID)
LarvByJarDat<- subset(LarvByJar, Block == "B1")
LarvByJarDat$JarID<- as.factor(LarvByJarDat$JarID)
LarvByJarDat$AccTankID<- as.factor(LarvByJarDat$AccTankID)
LarvByJarDat$TrtTankID<- as.factor(LarvByJarDat$TrtTankID)
LarvByJarDat$JarSeatable<- as.factor(LarvByJarDat$JarSeatable)
LarvByJarDat$ParentTrt<- as.factor(LarvByJarDat$ParentTrt)
LarvByJarDat$FemaleID<- as.factor(LarvByJarDat$FemaleID)
LarvByJarDat$MaleID<- as.factor(LarvByJarDat$MaleID)
LarvByJarDat$CrossID<- as.factor(LarvByJarDat$CrossID)
LarvByJarDat$JarTrt<- as.factor(LarvByJarDat$JarTrt)
#make a column for the ratio of major to minor axis
LarvaeDat$MajMinRat<- LarvaeDat$LarvaeMajorAxisLengthPix/LarvaeDat$LarvaeMinorAxisLengthPix
LarvaeDat$PerimDiamRat<-LarvaeDat$LarvaePerimeterum/LarvaeDat$LarvaeDiamumAnalyze jar chemistry data
CalJar<- subset(LarvByJarDat, JarOmegaCalcite >0)
CalJar<- subset(CalJar, BottleChemistry=="TRUE" )#only use those that had their chem measured
#look at the saturation state for jars
boxplot(JarOmegaCalcite~JarTrt*JarSeatable, data=CalJar)
#I was having trouble with making this model: JarOmegaCalcite~JarTrt + (1|ParentTrt)+ (1|TimeFilter), so I looked at TimeFilter because I think that is the problem.
Jar1<- lm(JarOmegaCalcite~TimeFilter, data=CalJar)
par(mfcol=c(2,2))
plot(Jar1)## Warning: not plotting observations with leverage one:
## 1, 2, 3, 4, 5, 6, 7, 8## Error in text.default(x, y, labels.id[ind], cex = cex.id, xpd = TRUE, : invalid 'pos' valuepar(mfcol=c(1,1))
#looks like we can't use timefilter because it is too specific to each jar. Get the error message "observations with leverage one"; try using seatable instead
Jar2<- lmer(JarOmegaCalcite~JarTrt + (1|JarSeatable), data=CalJar)
Jar2step<- step(Jar2)## Error in as_lmerModLmerTest(model): model not of class 'lmerMod': cannot coerce to class 'lmerModLmerTest#get error model not of class lmerMod. I looked up this issue and turns out that error results when all random effects are dropped. See https://stackoverflow.com/questions/55638476/how-to-do-stepwise-model-with-random-effect-lme4-lmertest for more info if you want.
JarFin<- lm(JarOmegaCalcite~JarTrt, data=CalJar)
#check assumptions
par(mfcol=c(2,2))
plot(JarFin)
par(mfcol=c(1,1))
acf(resid(JarFin))
summary(JarFin)##
## Call:
## lm(formula = JarOmegaCalcite ~ JarTrt, data = CalJar)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.11383 -0.02406 -0.01337 0.02975 0.19438
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.0863 0.0413 26.30 4.69e-09 ***
## JarTrtExposed -0.8486 0.0584 -14.53 4.93e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.09234 on 8 degrees of freedom
## Multiple R-squared: 0.9635, Adjusted R-squared: 0.9589
## F-statistic: 211.1 on 1 and 8 DF, p-value: 4.933e-07#Conclusion: JarSeatable doesn't need to be accounted for as a random effect in terms of water chem. Test for significant effect of parent treatment on egg size. Start with measured water chem
CalEggs<- subset(EggsB1, WCa_pHDIC != "NA")
CalEggs$TrtTankID<-droplevels(CalEggs)$TrtTankID
#start with models that have the response variable EggDiamum
eggdiam1<- lmer(EggDiamum~WCa_pHDIC*AdultLength + (1|FemaleID) + (1|TrtTankID) + (1|AccTankID), data=CalEggs)## boundary (singular) fit: see ?isSingulareggstep<-step(eggdiam1)## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingularprint(eggstep)## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 8 -1210.6 2437.2
## (1 | TrtTankID) 1 7 -1210.6 2435.2 0.000 1 1
## (1 | AccTankID) 2 6 -1210.6 2433.2 0.000 1 1
## (1 | FemaleID) 0 5 -1236.3 2482.7 51.425 1 7.439e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## WCa_pHDIC:AdultLength 1 3.5785 3.5785 1 11.920 0.1926 0.6686
## AdultLength 2 9.2137 9.2137 1 12.932 0.4958 0.4938
## WCa_pHDIC 3 23.0328 23.0328 1 14.281 1.2394 0.2840
##
## Model found:
## EggDiamum ~ (1 | FemaleID)#final model chosen: EggDiamum~ 1 + (1|FemaleID)
eggdiamfin<- lmer(EggDiamum~1+(1|FemaleID), data=CalEggs)
#Check assumptions
par(mfcol=c(1,1))
qqnorm(resid(eggdiamfin))
qqline(resid(eggdiamfin))
plot(eggdiamfin)
acf(resid(eggdiamfin))
#the above looks fine, seems to meet all assumptions.
summary(eggdiamfin)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: EggDiamum ~ 1 + (1 | FemaleID)
## Data: CalEggs
##
## REML criterion at convergence: 2431.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.8207 -0.6329 -0.1567 0.5843 3.4110
##
## Random effects:
## Groups Name Variance Std.Dev.
## FemaleID (Intercept) 4.331 2.081
## Residual 18.586 4.311
## Number of obs: 417, groups: FemaleID, 16
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 62.8803 0.5638 15.1316 111.5 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Conclusion: Female matters but not treatment or adult length. Acclimation tank and treatment tank also don't matter.
plot(EggDiamum~FemaleID, data=CalEggs) #we need to account for FemaleID, but not egg size in Larvae models. Can also make parallel models one with EggDiamum as a covariate and one with Female ID as a random effect
eggdiam1f<- lmer(EggDiamum~ParentTrt*AdultLength + (1|FemaleID) + (1|TrtTankID)+ (1|AccTankID), data=CalEggs)## boundary (singular) fit: see ?isSingulareggstepf<-step(eggdiam1f)## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingularprint(eggstepf)## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 8 -1211.0 2438.1
## (1 | TrtTankID) 1 7 -1211.0 2436.1 0.000 1 1
## (1 | AccTankID) 2 6 -1211.0 2434.1 0.000 1 1
## (1 | FemaleID) 0 5 -1236.9 2483.8 51.771 1 6.236e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## ParentTrt:AdultLength 1 4.6817 4.6817 1 11.912 0.2519 0.6249
## AdultLength 2 9.0581 9.0581 1 12.927 0.4874 0.4974
## ParentTrt 3 20.2635 20.2635 1 14.231 1.0904 0.3138
##
## Model found:
## EggDiamum ~ (1 | FemaleID)#final model chosen: EggDiamum~ 1 + (1|FemaleID), which is the same as the eggdiamfin model
summary(eggdiamfin)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: EggDiamum ~ 1 + (1 | FemaleID)
## Data: CalEggs
##
## REML criterion at convergence: 2431.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.8207 -0.6329 -0.1567 0.5843 3.4110
##
## Random effects:
## Groups Name Variance Std.Dev.
## FemaleID (Intercept) 4.331 2.081
## Residual 18.586 4.311
## Number of obs: 417, groups: FemaleID, 16
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 62.8803 0.5638 15.1316 111.5 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Choosing mean or median for larvae length data
ggplot(LarvaeDat, aes(x=LarvaeDiamum, color=JarTrt))+
geom_histogram(alpha=0.5, position="identity") +
facet_grid(~ParentTrt)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
larvmeans<- aggregate(LarvaeDiamum~ParentTrt*JarTrt, data=LarvaeDat, FUN=mean)#CP/CL : 75.24678; EP/CL: 75.46671; CP/EL: 64.22775; EP/EL: 65.62240
larvmedians<- aggregate(LarvaeDiamum~ParentTrt*JarTrt, data=LarvaeDat, FUN=median)#CP/CL : 75.29579; EP/CL: 75.52886; CP/EL: 64.27797; EP/EL: 65.82983
view(larvmeans)
view(larvmedians)
#I think we can use the mean for these data Look at the larvae diameter data
boxplot(LarvaeDiamum~FemaleID, data=LarvaeDat) #I looked specifically for weird things with EF07_B1 because that individual only had 10 eggs after removing outliers. But all seems fine! 
#create most complex model and then use the step function to see which model is best fit for the data; these are for ParentTrt and JarTrt as covariates
larvlen<-lmer(LarvaeDiamum~WCa_pHDIC*JarOmegaCalcite+ (1|FemaleID)+(1|MaleID)+(1|CrossID)+(1|JarID)+(1|JarSeatable)+(1|AccTankID), data=LarvaeDat)
steppedlarvlen<-step(larvlen)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00444468 (tol = 0.002, component 1)print(steppedlarvlen)## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 11 -6759.2 13540
## (1 | AccTankID) 1 10 -6759.5 13539 0.698 1 0.4034678
## (1 | CrossID) 2 9 -6760.6 13539 2.065 1 0.1507047
## (1 | FemaleID) 0 8 -6777.8 13572 34.482 1 4.301e-09 ***
## (1 | MaleID) 0 8 -6766.8 13550 12.376 1 0.0004349 ***
## (1 | JarID) 0 8 -6856.8 13730 192.380 1 < 2.2e-16 ***
## (1 | JarSeatable) 0 8 -6763.7 13543 6.191 1 0.0128429 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value
## WCa_pHDIC:JarOmegaCalcite 0 86.778 86.778 1 240.67 11.23
## Pr(>F)
## WCa_pHDIC:JarOmegaCalcite 0.0009345 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## LarvaeDiamum ~ WCa_pHDIC + JarOmegaCalcite + (1 | FemaleID) +
## (1 | MaleID) + (1 | JarID) + (1 | JarSeatable) + WCa_pHDIC:JarOmegaCalcite#final model chosen by the lme4 step function: y~WCa_pHDIC*JarOmegaCalcite +(1|FemaleID)+(1|MaleID)+(1|JarID)+ (1|JarSeatable)
larvlenfin<- lmer(LarvaeDiamum~WCa_pHDIC*JarOmegaCalcite+ (1|FemaleID)+(1|MaleID)+(1|JarID)+ (1|JarSeatable), data=LarvaeDat)
#check the final model to see that it meets assumptions
plot(larvlenfin) #seems to meet assumption of linearity and homoscedasticity
qqnorm(resid(larvlenfin))
qqline(resid(larvlenfin)) #has relatively long tails, doesn't seem to meet assumption of normality. Check on this
acf(resid(larvlenfin)) #this looks good 
summary(larvlenfin)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: LarvaeDiamum ~ WCa_pHDIC * JarOmegaCalcite + (1 | FemaleID) +
## (1 | MaleID) + (1 | JarID) + (1 | JarSeatable)
## Data: LarvaeDat
##
## REML criterion at convergence: 13521.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.8987 -0.5441 0.0495 0.6065 4.8949
##
## Random effects:
## Groups Name Variance Std.Dev.
## JarID (Intercept) 1.6344 1.2784
## FemaleID (Intercept) 0.7649 0.8746
## MaleID (Intercept) 0.3448 0.5872
## JarSeatable (Intercept) 0.1350 0.3675
## Residual 7.7275 2.7798
## Number of obs: 2698, groups:
## JarID, 270; FemaleID, 15; MaleID, 11; JarSeatable, 3
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 63.7296 0.7644 26.2756 83.369 < 2e-16 ***
## WCa_pHDIC -2.1627 0.8158 24.3372 -2.651 0.013893 *
## JarOmegaCalcite 10.8411 0.4941 240.7074 21.940 < 2e-16 ***
## WCa_pHDIC:JarOmegaCalcite 1.8102 0.5402 240.6723 3.351 0.000935 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) WC_pHDIC JrOmgC
## WCa_pHDIC -0.854
## JarOmegClct -0.428 0.391
## WC_HDIC:JOC 0.381 -0.437 -0.893#Conclusion: signficant main effect of Parent and Jar calcite sat, and sig. interaction between parent and Jar saturation statesUse fixed effect for Jar treatment and adult tank as covariate
larvlenJF<-lmer(LarvaeDiamum~WCa_pHDIC*JarTrt+ (1|FemaleID)+(1|MaleID)+(1|CrossID)+(1|JarID)+(1|JarSeatable)+(1|AccTankID), data=LarvaeDat)
steppedlarvlenJF<-step(larvlenJF)
print(steppedlarvlenJF)## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 11 -6757.4 13537
## (1 | AccTankID) 1 10 -6757.7 13535 0.702 1 0.402062
## (1 | CrossID) 2 9 -6758.5 13535 1.634 1 0.201178
## (1 | FemaleID) 0 8 -6776.1 13568 35.051 1 3.212e-09 ***
## (1 | MaleID) 0 8 -6766.1 13548 15.223 1 9.555e-05 ***
## (1 | JarID) 0 8 -6850.2 13716 183.331 1 < 2.2e-16 ***
## (1 | JarSeatable) 0 8 -6762.7 13541 8.292 1 0.003982 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## WCa_pHDIC:JarTrt 0 84.857 84.857 1 240.75 10.981 0.001062 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## LarvaeDiamum ~ WCa_pHDIC + JarTrt + (1 | FemaleID) + (1 | MaleID) +
## (1 | JarID) + (1 | JarSeatable) + WCa_pHDIC:JarTrt#final model chosen by the lme4 step function: y~WCa_pHDIC*JarTrt +(1|FemaleID)+(1|MaleID)+(1|JarID)+ (1|JarSeatable)
larvlenfinJF<- lmer(LarvaeDiamum~WCa_pHDIC*JarTrt+ (1|FemaleID)+(1|MaleID)+(1|JarID)+ (1|JarSeatable), data=LarvaeDat)
#check the final model to see that it meets assumptions
plot(larvlenfinJF) #seems to meet assumption of linearity and homoscedasticity
qqnorm(resid(larvlenfinJF))
qqline(resid(larvlenfinJF)) #has relatively long tails, I think this should be acceptable
acf(resid(larvlenfinJF)) #this looks good 
summary(larvlenfinJF)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: LarvaeDiamum ~ WCa_pHDIC * JarTrt + (1 | FemaleID) + (1 | MaleID) +
## (1 | JarID) + (1 | JarSeatable)
## Data: LarvaeDat
##
## REML criterion at convergence: 13517.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.9046 -0.5422 0.0507 0.6106 4.9021
##
## Random effects:
## Groups Name Variance Std.Dev.
## JarID (Intercept) 1.5778 1.2561
## FemaleID (Intercept) 0.7587 0.8711
## MaleID (Intercept) 0.3920 0.6261
## JarSeatable (Intercept) 0.1665 0.4081
## Residual 7.7275 2.7798
## Number of obs: 2698, groups:
## JarID, 270; FemaleID, 15; MaleID, 11; JarSeatable, 3
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 75.5290 0.7392 20.9914 102.172 < 2e-16 ***
## WCa_pHDIC -0.2155 0.7824 19.3222 -0.275 0.78593
## JarTrtExposed -9.2409 0.4149 240.7741 -22.273 < 2e-16 ***
## WCa_pHDIC:JarTrtExposed -1.5027 0.4535 240.7464 -3.314 0.00106 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) WC_pHDIC JrTrtE
## WCa_pHDIC -0.841
## JarTrtExpsd -0.281 0.259
## WC_HDIC:JTE 0.251 -0.292 -0.893Use fixed effect for both adult and larvae treatments
larvlenFE<-lmer(LarvaeDiamum~ParentTrt*JarTrt+ (1|FemaleID)+(1|MaleID)+(1|CrossID)+(1|JarID)+(1|JarSeatable)+(1|TrtTankID)+(1|AccTankID), data=LarvaeDat)## boundary (singular) fit: see ?isSingularsteppedlarvlenFE<-step(larvlenFE)## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00645175 (tol = 0.002, component 1)## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingularprint(steppedlarvlenFE)## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 12 -6758.3 13541
## (1 | TrtTankID) 1 11 -6758.3 13539 0.000 1 0.9995923
## (1 | AccTankID) 2 10 -6758.7 13537 0.641 1 0.4232503
## (1 | CrossID) 3 9 -6759.4 13537 1.580 1 0.2087381
## (1 | FemaleID) 0 8 -6776.6 13569 34.286 1 4.758e-09 ***
## (1 | MaleID) 0 8 -6767.0 13550 15.075 1 0.0001033 ***
## (1 | JarID) 0 8 -6852.0 13720 185.190 1 < 2.2e-16 ***
## (1 | JarSeatable) 0 8 -6763.6 13543 8.235 1 0.0041095 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## ParentTrt:JarTrt 0 75.391 75.391 1 240.74 9.7562 0.002006 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## LarvaeDiamum ~ ParentTrt + JarTrt + (1 | FemaleID) + (1 | MaleID) +
## (1 | JarID) + (1 | JarSeatable) + ParentTrt:JarTrt#final model chosen by the lme4 step function: y~ParentTrt*JarTrt +(1|FemaleID)+(1|MaleID)+(1|JarID)+ (1|JarSeatable)
larvlenfinFE<- lmer(LarvaeDiamum~ParentTrt*JarTrt+ (1|FemaleID)+(1|MaleID)+(1|JarID)+ (1|JarSeatable), data=LarvaeDat)
#check the final model to see that it meets assumptions
plot(larvlenfinFE) #seems to meet assumption of linearity and homoscedasticity
qqnorm(resid(larvlenfinFE))
qqline(resid(larvlenfinFE)) #has relatively long tails, I think this should be acceptable
acf(resid(larvlenfinFE)) #this looks good 
summary(larvlenfinFE)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: LarvaeDiamum ~ ParentTrt * JarTrt + (1 | FemaleID) + (1 | MaleID) +
## (1 | JarID) + (1 | JarSeatable)
## Data: LarvaeDat
##
## REML criterion at convergence: 13518.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.9019 -0.5407 0.0506 0.6101 4.9087
##
## Random effects:
## Groups Name Variance Std.Dev.
## JarID (Intercept) 1.5895 1.2607
## FemaleID (Intercept) 0.7497 0.8659
## MaleID (Intercept) 0.3909 0.6252
## JarSeatable (Intercept) 0.1663 0.4078
## Residual 7.7275 2.7798
## Number of obs: 2698, groups:
## JarID, 270; FemaleID, 15; MaleID, 11; JarSeatable, 3
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 75.2408 0.5109 15.6742 147.263 < 2e-16 ***
## ParentTrt2800 0.2381 0.6448 19.2484 0.369 0.71593
## JarTrtExposed -11.0160 0.2564 240.7031 -42.959 < 2e-16 ***
## ParentTrt2800:JarTrtExposed 1.1720 0.3752 240.7394 3.123 0.00201 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) PrT2800 JrTrtE
## PrntTrt2800 -0.624
## JarTrtExpsd -0.254 0.202
## PrT2800:JTE 0.174 -0.293 -0.683Visualize the larvae length data
#plot the data using the means for each jar
meanlarvjar<- ddply(LarvaeDat, .(JarID), numcolwise(mean, na.rm=T))
selarvjar<- ddply(LarvaeDat, .(JarID), numcolwise(se, na.rm=T))
ggplot(meanlarvjar,
aes(fill = WCa_pHDIC, y = LarvaeDiamum, x = JarOmegaCalcite)) +
geom_errorbar(aes(ymin= meanlarvjar$LarvaeDiamum-selarvjar$LarvaeDiamum, ymax=meanlarvjar$LarvaeDiamum+selarvjar$LarvaeDiamum),width=0.05, position=position_dodge(0))+
geom_point(aes(colour=WCa_pHDIC)) +
ggtitle("Larvae diameter by larvae jar saturation state")+
ylab("Larvae Diameter (um)")+
xlab("Larvae Jar Calcite Saturation State")+
theme_classic()## Warning: Use of `meanlarvjar$LarvaeDiamum` is discouraged. Use `LarvaeDiamum`
## instead.
## Warning: Use of `meanlarvjar$LarvaeDiamum` is discouraged. Use `LarvaeDiamum`
## instead.## Warning: position_dodge requires non-overlapping x intervals
#the reason there is spread is that for the jars that we actually measured carbonate dynamics I used those values, for all other jars, I used the average of those measurements
#Try to visualize it in a more pleasant way
#get means and SE of larvae diameters
meanlarv<- ddply(LarvaeDat, .(JarTrt, ParentTrt), numcolwise(mean, na.rm=T))
selarv<- ddply(LarvaeDat, .(JarTrt, ParentTrt), numcolwise(se, na.rm=T))
#plot the mean Larvae Diameter +/- SE
ggplot(data = meanlarv, aes(x = JarOmegaCalcite, y = LarvaeDiamum, group=ParentTrt, fill=ParentTrt)) +
geom_errorbar(aes(ymin= LarvaeDiamum-selarv$LarvaeDiamum, ymax=LarvaeDiamum+selarv$LarvaeDiamum),width=0.05, position=position_dodge(0))+
geom_errorbarh(aes(xmin= JarOmegaCalcite-selarv$JarOmegaCalcite, xmax=JarOmegaCalcite+selarv$JarOmegaCalcite),width=0.05, position=position_dodge(0))+
geom_point(aes(colour=ParentTrt))+
labs(x = "Larvae Jar Saturation State (Calcite)", y = "Larvae Diameter (um)") +
scale_color_manual(values = c("skyblue2", "red2"), name="Parent Treatment") +
theme_classic()## Warning: Ignoring unknown parameters: width
## Warning: position_dodge requires non-overlapping x intervals
## Warning: position_dodge requires non-overlapping x intervals
#alternatively, I could visualize with bargraphs:
#plot with parent treatment on the x-axis
ggplot(data = meanlarv, aes(x = ParentTrt , y = LarvaeDiamum, fill=JarTrt)) +
geom_bar(stat="identity", aes(fill=JarTrt), position="dodge")+
labs(x = "Parent Treatment", y = "Larvae Diameter (um)") +
scale_fill_manual(name = "Larvae Treatment", labels = c("Control", "OA Exposed"), values = c("skyblue2", "red2")) +
geom_errorbar(aes(ymin= meanlarv$LarvaeDiamum-selarv$LarvaeDiamum, ymax=meanlarv$LarvaeDiamum+selarv$LarvaeDiamum),width=0.2, position=position_dodge(.9))+
theme_classic()## Warning: Use of `meanlarv$LarvaeDiamum` is discouraged. Use `LarvaeDiamum`
## instead.## Warning: Use of `meanlarv$LarvaeDiamum` is discouraged. Use `LarvaeDiamum`
## instead.
#now plot the difference between parental treatments for larvae diameter
crossdiff <- data.frame(tapply(LarvaeDat$LarvaeDiamum, list(LarvaeDat$CrossID, LarvaeDat$JarTrt), mean)) %>%
rownames_to_column("CrossID") %>%
mutate(CrossTemp = CrossID) %>%
mutate(ParentTrt = ifelse(startsWith(CrossID, "E"), "Exposed", "Control")) %>%
separate(CrossTemp, c("FemaleID", "MaleID", "BlockID", "Fourth"), sep="_") %>%
mutate(Diff = Exposed - Control)
crossdiff$CrossID<- as.factor(crossdiff$CrossID)
crossdiff$FemaleID<-as.factor(crossdiff$FemaleID)
crossdiff$MaleID<-as.factor(crossdiff$MaleID)
#look at cross differences by Male and Female ID
ggplot(data = crossdiff, aes(x = MaleID, y = Diff)) +
geom_point(aes(colour=FemaleID)) +
labs(x = "Male ID", y = "Change in mean shell length (um) per family\nfrom control to OA conditions") +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), panel.background = element_blank(),
axis.line = element_line(colour = "black"), panel.border=element_rect(colour="black", fill=NA))+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5))
#get the means of the differences to visualize overall differences
meandiff<- aggregate(Diff~ParentTrt, data=crossdiff, FUN=mean)
sediff<- aggregate(Diff~ParentTrt, data=crossdiff, FUN=se)
ggplot(data = meandiff, aes(x = ParentTrt, y = Diff)) +
geom_bar(stat="identity", aes(fill="gray45"), position="dodge") +
labs(x = "Parental Treatment", y = "Change in mean shell length (um) per family\nfrom control to OA conditions") +
scale_x_discrete(breaks = c("Control", "Exposed"), labels = c("Control", "OA Exposed")) +
scale_fill_manual(values = "gray45") +
geom_errorbar(aes(ymin= Diff-sediff$Diff, ymax=Diff+sediff$Diff),width=0.2, position=position_dodge(.9))+
theme(legend.position = "none", panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), panel.background = element_blank(),
axis.line = element_line(colour = "black"))
diffttest<- t.test(Diff~ParentTrt, data=crossdiff, var.equal=T)
diffttest #mean shell length difference is significantly higher for OA parent treatment##
## Two Sample t-test
##
## data: Diff by ParentTrt
## t = -2.6006, df = 43, p-value = 0.01271
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -2.1378060 -0.2703394
## sample estimates:
## mean in group Control mean in group Exposed
## -11.049847 -9.845775Analyze at growth per day. This will likely be the main figure
boxplot(GrowthPerDay~ParentTrt*JarTrt, data=LarvaeDat)
larvgrow<-lmer(GrowthPerDay~WCa_pHDIC*JarOmegaCalcite+ (1|FemaleID)+(1|MaleID)+(1|CrossID)+(1|JarID) +(1|JarSeatable)+(1|AccTankID) , data=LarvaeDat)## boundary (singular) fit: see ?isSingularsteppedlarvgrow<- step(larvgrow)## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## unable to evaluate scaled gradient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge: degenerate Hessian with 1 negative eigenvalues## Warning: Model failed to converge with 1 negative eigenvalue: -1.8e+00print(steppedlarvgrow)## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 11 -3809.8 7641.7
## (1 | AccTankID) 1 10 -3809.8 7639.7 0.000 1 0.999686
## (1 | FemaleID) 0 9 -3828.6 7675.1 37.443 1 9.411e-10 ***
## (1 | MaleID) 0 9 -3811.7 7641.5 3.782 1 0.051815 .
## (1 | CrossID) 0 9 -3814.3 7646.7 9.005 1 0.002692 **
## (1 | JarID) 0 9 -3891.5 7801.0 163.351 1 < 2.2e-16 ***
## (1 | JarSeatable) 0 9 -3813.1 7644.1 6.449 1 0.011100 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value
## WCa_pHDIC:JarOmegaCalcite 0 10.242 10.242 1 221.49 11.928
## Pr(>F)
## WCa_pHDIC:JarOmegaCalcite 0.0006624 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## GrowthPerDay ~ WCa_pHDIC + JarOmegaCalcite + (1 | FemaleID) +
## (1 | MaleID) + (1 | CrossID) + (1 | JarID) + (1 | JarSeatable) +
## WCa_pHDIC:JarOmegaCalcite#final model chosen by the lme4 step function: y~WCa_pHDIC*JarOmegaCalcite+(1|FemaleID)+(1|MaleID)+ (1|CrossID)+ (1|JarID) + (1|JarSeatable)
larvgrowfin<- lmer(GrowthPerDay~WCa_pHDIC*JarOmegaCalcite+(1|FemaleID)+(1|MaleID)+(1|CrossID) + (1|JarID) + (1|JarSeatable), data=LarvaeDat)
#check the final model to see that it meets assumptions
plot(larvgrowfin) #seems to meet assumption of linearity and homoscedasticity
qqnorm(resid(larvgrowfin))
qqline(resid(larvgrowfin)) #has relatively long tails, based on scale I think it reasonably meets assumption of normality. 
acf(resid(larvgrowfin))
summary(larvgrowfin)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: GrowthPerDay ~ WCa_pHDIC * JarOmegaCalcite + (1 | FemaleID) +
## (1 | MaleID) + (1 | CrossID) + (1 | JarID) + (1 | JarSeatable)
## Data: LarvaeDat
##
## REML criterion at convergence: 7619.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.9566 -0.5404 0.0513 0.6056 4.9473
##
## Random effects:
## Groups Name Variance Std.Dev.
## JarID (Intercept) 0.16985 0.4121
## CrossID (Intercept) 0.02121 0.1456
## FemaleID (Intercept) 0.43849 0.6622
## MaleID (Intercept) 0.03583 0.1893
## JarSeatable (Intercept) 0.01485 0.1218
## Residual 0.85862 0.9266
## Number of obs: 2698, groups:
## JarID, 270; CrossID, 45; FemaleID, 15; MaleID, 11; JarSeatable, 3
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.8459 0.4266 17.7762 1.983 0.063027 .
## WCa_pHDIC -1.3679 0.4621 17.1680 -2.960 0.008705 **
## JarOmegaCalcite 3.6130 0.1611 221.3094 22.432 < 2e-16 ***
## WCa_pHDIC:JarOmegaCalcite 0.6082 0.1761 221.4933 3.454 0.000662 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) WC_pHDIC JrOmgC
## WCa_pHDIC -0.879
## JarOmegClct -0.250 0.225
## WC_HDIC:JOC 0.223 -0.251 -0.893#Conclusion: Signficant main effect of tank and jar calcite saturation state and significant interaction.Analyze larvae growth with adult trt as covariate and jar as fixed effect
larvgrowJF<-lmer(GrowthPerDay~WCa_pHDIC*JarTrt+ (1|FemaleID)+(1|MaleID)+(1|CrossID)+(1|JarID) +(1|JarSeatable)+(1|AccTankID) , data=LarvaeDat)## boundary (singular) fit: see ?isSingularsteppedlarvgrowJF<- step(larvgrowJF)## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## unable to evaluate scaled gradient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge: degenerate Hessian with 1 negative eigenvalues## Warning: Model failed to converge with 1 negative eigenvalue: -1.7e+00print(steppedlarvgrowJF)## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 11 -3808.0 7638.1
## (1 | AccTankID) 1 10 -3808.0 7636.1 0.000 1 1.000000
## (1 | FemaleID) 0 9 -3827.0 7672.0 37.986 1 7.126e-10 ***
## (1 | MaleID) 0 9 -3810.6 7639.1 5.058 1 0.024518 *
## (1 | CrossID) 0 9 -3812.3 7642.5 8.449 1 0.003653 **
## (1 | JarID) 0 9 -3886.6 7791.2 157.122 1 < 2.2e-16 ***
## (1 | JarSeatable) 0 9 -3812.3 7642.6 8.566 1 0.003426 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## WCa_pHDIC:JarTrt 0 9.9793 9.9793 1 221.43 11.623 0.0007741 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## GrowthPerDay ~ WCa_pHDIC + JarTrt + (1 | FemaleID) + (1 | MaleID) +
## (1 | CrossID) + (1 | JarID) + (1 | JarSeatable) + WCa_pHDIC:JarTrt#final model chosen by the lme4 step function: y~WCa_pHDIC*JarTrt+(1|FemaleID)+(1|MaleID)+ (1|CrossID)+ (1|JarID) + (1|JarSeatable)
larvgrowfinJF<- lmer(GrowthPerDay~WCa_pHDIC*JarTrt+(1|FemaleID)+(1|MaleID)+(1|CrossID) + (1|JarID) + (1|JarSeatable), data=LarvaeDat)
#check the final model to see that it meets assumptions
plot(larvgrowfinJF) #seems to meet assumption of linearity and homoscedasticity
qqnorm(resid(larvgrowfinJF))
qqline(resid(larvgrowfinJF)) #has relatively long tails, based on scale I think it reasonably meets assumption of normality. 
acf(resid(larvgrowfinJF))
summary(larvgrowfinJF)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: GrowthPerDay ~ WCa_pHDIC * JarTrt + (1 | FemaleID) + (1 | MaleID) +
## (1 | CrossID) + (1 | JarID) + (1 | JarSeatable)
## Data: LarvaeDat
##
## REML criterion at convergence: 7616.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.9576 -0.5394 0.0530 0.6100 4.9506
##
## Random effects:
## Groups Name Variance Std.Dev.
## JarID (Intercept) 0.16522 0.4065
## CrossID (Intercept) 0.01806 0.1344
## FemaleID (Intercept) 0.43793 0.6618
## MaleID (Intercept) 0.04226 0.2056
## JarSeatable (Intercept) 0.01830 0.1353
## Residual 0.85862 0.9266
## Number of obs: 2698, groups:
## JarID, 270; CrossID, 45; FemaleID, 15; MaleID, 11; JarSeatable, 3
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 4.7790 0.4222 17.0367 11.319 2.39e-09 ***
## WCa_pHDIC -0.7149 0.4565 16.3546 -1.566 0.136483
## JarTrtExposed -3.0787 0.1356 221.2118 -22.703 < 2e-16 ***
## WCa_pHDIC:JarTrtExposed -0.5054 0.1482 221.4272 -3.409 0.000774 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) WC_pHDIC JrTrtE
## WCa_pHDIC -0.876
## JarTrtExpsd -0.161 0.145
## WC_HDIC:JTE 0.144 -0.163 -0.893Analyze larvae growth with adult trt as covariate and jar as fixed effect
larvgrowFE<-lmer(GrowthPerDay~ParentTrt*JarTrt+ (1|FemaleID)+(1|MaleID)+(1|CrossID)+(1|JarID) +(1|JarSeatable)+(1|TrtTankID)+(1|AccTankID) , data=LarvaeDat)## boundary (singular) fit: see ?isSingularsteppedlarvgrowFE<- step(larvgrowFE)## boundary (singular) fit: see ?isSingular## Warning: Model failed to converge with 1 negative eigenvalue: -9.5e+01## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular## Warning: Model failed to converge with 1 negative eigenvalue: -9.9e+00## boundary (singular) fit: see ?isSingular## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.113534 (tol = 0.002, component 1)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00209849 (tol = 0.002, component 1)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.113534 (tol = 0.002, component 1)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00235773 (tol = 0.002, component 1)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.006054 (tol = 0.002, component 1)print(steppedlarvgrowFE)## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 12 -3809.2 7642.4
## (1 | AccTankID) 1 11 -3809.2 7640.4 0.000 1 1.000000
## (1 | TrtTankID) 2 10 -3809.2 7638.4 0.000 1 0.999366
## (1 | CrossID) 3 9 -3809.9 7637.7 1.310 1 0.252349
## (1 | FemaleID) 0 8 -3905.5 7826.9 191.164 1 < 2.2e-16 ***
## (1 | MaleID) 0 8 -3818.2 7652.4 16.666 1 4.458e-05 ***
## (1 | JarID) 0 8 -3902.2 7820.4 184.710 1 < 2.2e-16 ***
## (1 | JarSeatable) 0 8 -3813.9 7643.9 8.156 1 0.004293 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## ParentTrt:JarTrt 0 8.3917 8.3917 1 241.78 9.7721 0.001989 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## GrowthPerDay ~ ParentTrt + JarTrt + (1 | FemaleID) + (1 | MaleID) +
## (1 | JarID) + (1 | JarSeatable) + ParentTrt:JarTrt#final model chosen by the lme4 step function: y~WCa_pHDIC*JarTrt+(1|FemaleID)+(1|MaleID)+ (1|CrossID)+ (1|JarID) + (1|JarSeatable)
larvgrowfinFE<- lmer(GrowthPerDay~ParentTrt*JarTrt+(1|FemaleID)+(1|MaleID)+(1|CrossID) + (1|JarID) + (1|JarSeatable), data=LarvaeDat)
#check the final model to see that it meets assumptions
plot(larvgrowfinFE) #seems to meet assumption of linearity and homoscedasticity
qqnorm(resid(larvgrowfinFE))
qqline(resid(larvgrowfinFE)) #has relatively long tails, based on scale I think it reasonably meets assumption of normality. 
acf(resid(larvgrowfinFE))
summary(larvgrowfinFE)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: GrowthPerDay ~ ParentTrt * JarTrt + (1 | FemaleID) + (1 | MaleID) +
## (1 | CrossID) + (1 | JarID) + (1 | JarSeatable)
## Data: LarvaeDat
##
## REML criterion at convergence: 7618.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.9541 -0.5394 0.0537 0.6073 4.9570
##
## Random effects:
## Groups Name Variance Std.Dev.
## JarID (Intercept) 0.16666 0.4082
## CrossID (Intercept) 0.01776 0.1333
## FemaleID (Intercept) 0.45139 0.6719
## MaleID (Intercept) 0.04228 0.2056
## JarSeatable (Intercept) 0.01828 0.1352
## Residual 0.85861 0.9266
## Number of obs: 2698, groups:
## JarID, 270; CrossID, 45; FemaleID, 15; MaleID, 11; JarSeatable, 3
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.92983 0.27447 18.11960 14.318 2.54e-11 ***
## ParentTrt2800 0.56798 0.38169 16.29530 1.488 0.15583
## JarTrtExposed -3.67574 0.08388 221.74326 -43.823 < 2e-16 ***
## ParentTrt2800:JarTrtExposed 0.39433 0.12271 221.44574 3.214 0.00151 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) PrT2800 JrTrtE
## PrntTrt2800 -0.661
## JarTrtExpsd -0.155 0.111
## PrT2800:JTE 0.106 -0.162 -0.684Larvae growth per day visualization
#plot growth per day
ggplot(LarvaeDat,
aes(fill = WCa_pHDIC, y = GrowthPerDay, x = JarOmegaCalcite)) +
geom_point(aes(colour=WCa_pHDIC)) +
ggtitle("Larvae growth by parent treatment")+
ylab("Growth (um/day)")+
theme_classic()
#get means and SE of larvae growth
meangrowth<- aggregate(GrowthPerDay~ParentTrt*JarTrt, data=LarvaeDat, FUN=mean)
segrowth<- aggregate(GrowthPerDay~ParentTrt*JarTrt, data=LarvaeDat, FUN=se)
meangrowth$JarTrt<- as.factor(meangrowth$JarTrt)
segrowth$JarTrt<- as.factor(segrowth$JarTrt)
#use a bargraph to plot the mean Larvae Diameter +/- SE
ggplot(data = meangrowth, aes(x = ParentTrt, y = GrowthPerDay, fill=JarTrt)) +
geom_bar(stat="identity", aes(fill=JarTrt), position="dodge")+
labs(x = "Parent Treatment", y = "Larvae Growth (um/day)") +
scale_fill_manual(name = "Larvae Treatment", labels = c("Control", "OA Exposed"), values = c("skyblue2", "red2")) +
geom_errorbar(aes(ymin= GrowthPerDay-segrowth$GrowthPerDay, ymax=GrowthPerDay+segrowth$GrowthPerDay),width=0.2, position=position_dodge(.9))+
theme_classic()
#make plot like the Chirgwin et al. paper
par(mar = c(5, 4, 0.5, 0.5), bg = "transparent")
barCenters <- barplot(meangrowth$GrowthPerDay~meangrowth$JarTrt+meangrowth$ParentTrt,
space= c(0,.25),
beside = TRUE, las = 1,
ylim = c(0, 5),
cex.names = 0.75,
ylab = expression(paste("Growth (", mu, "m/day)")),
xlab = "Larvae Treatment",
border = "black", axes = TRUE,
names=c("Control", "OA Exposed", "Control", "OA Exposed"),
legend.text = FALSE, col= c(rgb(0, 0, 1, alpha = 0.5), rgb(1, 0, 0, alpha = 0.5)))
arrows(x0=barCenters, y0=meangrowth$GrowthPerDay[c(1,3,2,4)] - segrowth$GrowthPerDay[c(1,3,2,4)], y1=meangrowth$GrowthPerDay[c(1,3,2,4)] + segrowth$GrowthPerDay[c(1,3,2,4)], angle=90, code=3)
rect(xleft = -0.2,
ybottom = 0,
xright = 2.37,
ytop = 5,
border = TRUE,
col = rgb(0, 0, 1, alpha = 0.15))
rect(xleft = 2.37,
ybottom = 0,
xright = 5.2,
ytop = 5,
border = TRUE,
col = rgb(1, 0, 0, alpha = 0.15))
barplot(meangrowth$GrowthPerDay~meangrowth$JarTrt+meangrowth$ParentTrt,
space= c(0,.25),
beside = TRUE, las = 1,
ylim = c(0, 5),
cex.names = 0.75,
border = "black", axes = FALSE,
names=c("Control", "OA Exposed", "Control", "OA Exposed"),
legend.text = FALSE, col= c(rgb(0, 0, 1, alpha = 0.5), rgb(1, 0, 0, alpha = 0.5)), add=TRUE)
#now plot the difference between parental treatments for larvae growth
growthdiff <- data.frame(tapply(LarvaeDat$GrowthPerDay, list(LarvaeDat$CrossID, LarvaeDat$JarTrt), mean)) %>%
rownames_to_column("CrossID") %>%
mutate(CrossTemp = CrossID) %>%
mutate(ParentTrt = ifelse(startsWith(CrossID, "E"), "Exposed", "Control")) %>%
separate(CrossTemp, c("FemaleID", "MaleID", "BlockID", "Fourth"), sep="_") %>% mutate(Diff = Exposed - Control)
growthdiff$CrossID<- as.factor(growthdiff$CrossID)
growthdiff$FemaleID<-as.factor(growthdiff$FemaleID)
growthdiff$MaleID<-as.factor(growthdiff$MaleID)
#look at cross differences by Male and Female ID
ggplot(data = growthdiff, aes(x = MaleID, y = Diff)) +
geom_point(aes(colour=FemaleID)) +
labs(x = "Male ID", y = "Change in larvae growth (um/day) per family\nfrom control to OA conditions") +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), panel.background = element_blank(),
axis.line = element_line(colour = "black"), panel.border=element_rect(colour="black", fill=NA))+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5))
#get the means of the differences to visualize overall differences
meangrdiff<- aggregate(Diff~ParentTrt, data=growthdiff, FUN=mean)
segrdiff<- aggregate(Diff~ParentTrt, data=growthdiff, FUN=se)
ggplot(data = meangrdiff, aes(x = ParentTrt, y = Diff)) +
geom_bar(stat="identity", aes(fill="gray45"), position="dodge") +
labs(x = "Parental Treatment", y = "Change in larvae growth (um) per family\nfrom control to OA conditions") +
scale_x_discrete(breaks = c("Control", "Exposed"), labels = c("Control", "OA")) +
scale_fill_manual(values = "gray45") +
geom_errorbar(aes(ymin= Diff-segrdiff$Diff, ymax=Diff+segrdiff$Diff),width=0.2, position=position_dodge(.9))+
theme(legend.position = "none", panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), panel.background = element_blank(),
axis.line = element_line(colour = "black"))
diffgrttest<- t.test(Diff~ParentTrt, data=growthdiff, var.equal=T)
diffgrttest #significant effect of parent treatment on growth. ##
## Two Sample t-test
##
## data: Diff by ParentTrt
## t = -2.6006, df = 43, p-value = 0.01271
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.71260201 -0.09011314
## sample estimates:
## mean in group Control mean in group Exposed
## -3.683282 -3.281925Look at larvae area
boxplot(LarvaeAreaum2~ParentTrt*JarTrt, data=LarvaeDat)
#create most complex model and then use the step function to see which model is best fit for the data; these are for ParentTrt and JarTrt as covariates
larvarea<-lmer(LarvaeAreaum2~WCa_pHDIC*JarOmegaCalcite+ (1|FemaleID)+(1|MaleID)+(1|CrossID)+(1|JarSeatable)+(1|JarID)+(1|AccTankID) , data=LarvaeDat)
steppedlarvarea<-step(larvarea)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00443797 (tol = 0.002, component 1)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00338512 (tol = 0.002, component 1)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00332707 (tol = 0.002, component 1)print(steppedlarvarea)## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 11 -18923 37868
## (1 | MaleID) 1 10 -18923 37867 1.097 1 0.2949993
## (1 | AccTankID) 2 9 -18924 37866 1.464 1 0.2262316
## (1 | JarSeatable) 3 8 -18925 37866 2.205 1 0.1375245
## (1 | FemaleID) 0 7 -18928 37870 5.405 1 0.0200845 *
## (1 | CrossID) 0 7 -18931 37876 12.013 1 0.0005284 ***
## (1 | JarID) 0 7 -19033 38079 214.833 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value
## WCa_pHDIC:JarOmegaCalcite 0 376864 376864 1 223.37 5.9083
## Pr(>F)
## WCa_pHDIC:JarOmegaCalcite 0.01586 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## LarvaeAreaum2 ~ WCa_pHDIC + JarOmegaCalcite + (1 | FemaleID) +
## (1 | CrossID) + (1 | JarID) + WCa_pHDIC:JarOmegaCalcite#final model chosen by the lme4 step function: y~WCa_pHDIC*JarOmegaCalcite+(1|FemaleID)+(1|CrossID)+(1|JarID)
larvareafin<- lmer(LarvaeAreaum2~WCa_pHDIC*JarOmegaCalcite+(1|FemaleID)+(1|CrossID)+(1|JarID), data=LarvaeDat)
#check the final model to see that it meets assumptions
plot(larvareafin) #seems to meet assumption of linearity and homoscedasticity
qqnorm(resid(larvareafin))
qqline(resid(larvareafin)) #has relatively long tails, doesn't seem to meet assumption of normality. 
acf(resid(larvareafin)) #this looks fine to me. 
summary(larvareafin) #try transforming the data## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: LarvaeAreaum2 ~ WCa_pHDIC * JarOmegaCalcite + (1 | FemaleID) +
## (1 | CrossID) + (1 | JarID)
## Data: LarvaeDat
##
## REML criterion at convergence: 37850.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.5267 -0.5798 0.0408 0.6257 5.2994
##
## Random effects:
## Groups Name Variance Std.Dev.
## JarID (Intercept) 15333 123.82
## CrossID (Intercept) 4909 70.06
## FemaleID (Intercept) 5197 72.09
## Residual 63786 252.56
## Number of obs: 2698, groups: JarID, 270; CrossID, 45; FemaleID, 15
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2660.72 60.11 24.02 44.267 <2e-16 ***
## WCa_pHDIC -151.83 65.62 23.94 -2.314 0.0296 *
## JarOmegaCalcite 1142.11 46.93 223.22 24.337 <2e-16 ***
## WCa_pHDIC:JarOmegaCalcite 124.73 51.32 223.37 2.431 0.0159 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) WC_pHDIC JrOmgC
## WCa_pHDIC -0.893
## JarOmegClct -0.517 0.462
## WC_HDIC:JOC 0.461 -0.516 -0.893larvareat<-lmer(sqrt(LarvaeAreaum2)~WCa_pHDIC*JarOmegaCalcite+ (1|FemaleID)+(1|MaleID)+(1|CrossID)+(1|JarSeatable)+(1|JarID)+(1|AccTankID) , data=LarvaeDat)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00224743 (tol = 0.002, component 1)steppedlarvareat<-step(larvareat)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0020834 (tol = 0.002, component 1)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00754904 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00754904 (tol = 0.002, component 1)print(steppedlarvareat)## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 11 -6182.9 12388
## (1 | AccTankID) 1 10 -6183.5 12387 1.118 1 0.2903515
## (1 | MaleID) 2 9 -6184.4 12387 1.825 1 0.1767765
## (1 | JarSeatable) 3 8 -6185.4 12387 2.114 1 0.1459162
## (1 | FemaleID) 0 7 -6188.1 12390 5.308 1 0.0212333 *
## (1 | CrossID) 0 7 -6191.3 12396 11.643 1 0.0006446 ***
## (1 | JarID) 0 7 -6286.0 12586 201.169 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value
## WCa_pHDIC:JarOmegaCalcite 0 35.903 35.903 1 223.38 7.1837
## Pr(>F)
## WCa_pHDIC:JarOmegaCalcite 0.007905 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## sqrt(LarvaeAreaum2) ~ WCa_pHDIC + JarOmegaCalcite + (1 | FemaleID) +
## (1 | CrossID) + (1 | JarID) + WCa_pHDIC:JarOmegaCalcite#final model chosen by the lme4 step function: y~WCa_pHDIC*JarOmegaCalcite+(1|FemaleID)+(1|CrossID)+(1|JarID)
larvareafint<- lmer(sqrt(LarvaeAreaum2)~WCa_pHDIC*JarOmegaCalcite+(1|FemaleID)+(1|CrossID)+(1|JarID), data=LarvaeDat)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00754904 (tol = 0.002, component 1)#got warning message that model failed to converged. I read this stack overflow talking about checking to see if it is acceptable: https://stats.stackexchange.com/questions/97929/lmer-model-fails-to-converge
TestConverge <- with(larvareafint@optinfo$derivs,solve(Hessian,gradient))
max(abs(TestConverge))## [1] 0.0004658432#this value seems small enough that it's acceptable
plot(larvareafint) #seems to meet assumption of linearity and homoscedasticity
qqnorm(resid(larvareafint))
qqline(resid(larvareafint)) #this seems to have improved it
acf(resid(larvareafint)) #this looks fine to me. 
summary(larvareafint)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: sqrt(LarvaeAreaum2) ~ WCa_pHDIC * JarOmegaCalcite + (1 | FemaleID) +
## (1 | CrossID) + (1 | JarID)
## Data: LarvaeDat
##
## REML criterion at convergence: 12370.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.5313 -0.5499 0.0578 0.6247 5.2230
##
## Random effects:
## Groups Name Variance Std.Dev.
## JarID (Intercept) 1.1442 1.0697
## CrossID (Intercept) 0.3637 0.6031
## FemaleID (Intercept) 0.3846 0.6202
## Residual 4.9979 2.2356
## Number of obs: 2698, groups: JarID, 270; CrossID, 45; FemaleID, 15
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 51.7442 0.5194 24.2032 99.624 < 2e-16 ***
## WCa_pHDIC -1.4302 0.5670 24.1239 -2.522 0.01866 *
## JarOmegaCalcite 9.8236 0.4084 223.2301 24.056 < 2e-16 ***
## WCa_pHDIC:JarOmegaCalcite 1.1968 0.4465 223.3805 2.680 0.00791 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) WC_pHDIC JrOmgC
## WCa_pHDIC -0.893
## JarOmegClct -0.520 0.465
## WC_HDIC:JOC 0.464 -0.519 -0.893
## convergence code: 0
## Model failed to converge with max|grad| = 0.00754904 (tol = 0.002, component 1)#significant main effects of tank and jar treatments and significant interactionAnalyze surface area with Parent Trt as covariate and Jar Treatment as fixed effect
larvareatJF<-lmer(sqrt(LarvaeAreaum2)~WCa_pHDIC*JarTrt+ (1|FemaleID)+(1|MaleID)+(1|CrossID)+(1|JarSeatable)+(1|JarID)+(1|AccTankID) , data=LarvaeDat)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00357282 (tol = 0.002, component 1)steppedlarvareatJF<-step(larvareatJF)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00296653 (tol = 0.002, component 1)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00461815 (tol = 0.002, component 1)print(steppedlarvareatJF)## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 11 -6181.1 12384
## (1 | AccTankID) 1 10 -6181.6 12383 1.111 1 0.2918610
## (1 | CrossID) 2 9 -6182.9 12384 2.640 1 0.1042159
## (1 | FemaleID) 0 8 -6201.7 12419 37.510 1 9.095e-10 ***
## (1 | MaleID) 0 8 -6190.0 12396 14.136 1 0.0001701 ***
## (1 | JarSeatable) 0 8 -6184.7 12385 3.444 1 0.0634889 .
## (1 | JarID) 0 8 -6292.4 12601 218.856 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## WCa_pHDIC:JarTrt 0 32.41 32.41 1 240.67 6.4848 0.0115 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## sqrt(LarvaeAreaum2) ~ WCa_pHDIC + JarTrt + (1 | FemaleID) + (1 |
## MaleID) + (1 | JarSeatable) + (1 | JarID) + WCa_pHDIC:JarTrt#final model chosen by the lme4 step function: y~WCa_pHDIC*JarTrt+(1|FemaleID)+(1|MaleID)+(1|JarID)+(1|JarSeatable)
larvareafintJF<- lmer(sqrt(LarvaeAreaum2)~WCa_pHDIC*JarTrt+(1|FemaleID)+(1|MaleID)+(1|JarID)+(1|JarSeatable), data=LarvaeDat)
plot(larvareafintJF) #seems to meet assumption of linearity and homoscedasticity
qqnorm(resid(larvareafintJF))
qqline(resid(larvareafintJF)) #this seems to have improved it
acf(resid(larvareafintJF)) #this looks fine to me. 
summary(larvareafintJF)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: sqrt(LarvaeAreaum2) ~ WCa_pHDIC * JarTrt + (1 | FemaleID) + (1 |
## MaleID) + (1 | JarID) + (1 | JarSeatable)
## Data: LarvaeDat
##
## REML criterion at convergence: 12365.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.4593 -0.5425 0.0603 0.6299 5.1752
##
## Random effects:
## Groups Name Variance Std.Dev.
## JarID (Intercept) 1.16250 1.0782
## FemaleID (Intercept) 0.55911 0.7477
## MaleID (Intercept) 0.26504 0.5148
## JarSeatable (Intercept) 0.06002 0.2450
## Residual 4.99788 2.2356
## Number of obs: 2698, groups:
## JarID, 270; FemaleID, 15; MaleID, 11; JarSeatable, 3
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 62.4349 0.6082 20.9581 102.663 <2e-16 ***
## WCa_pHDIC -0.1523 0.6598 19.3614 -0.231 0.8199
## JarTrtExposed -8.3685 0.3489 240.6878 -23.986 <2e-16 ***
## WCa_pHDIC:JarTrtExposed -0.9711 0.3813 240.6713 -2.547 0.0115 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) WC_pHDIC JrTrtE
## WCa_pHDIC -0.863
## JarTrtExpsd -0.287 0.259
## WC_HDIC:JTE 0.257 -0.291 -0.893Analyze surface area with Parent and Jar Treatment as fixed effects
larvareatFE<-lmer(sqrt(LarvaeAreaum2)~ParentTrt*JarTrt+ (1|FemaleID)+(1|MaleID)+(1|CrossID)+(1|JarSeatable)+(1|JarID)+(1|AccTankID) + (1|TrtTankID), data=LarvaeDat)## boundary (singular) fit: see ?isSingularsteppedlarvareatFE<-step(larvareatFE)## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00347465 (tol = 0.002, component 1)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00448054 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00448054 (tol = 0.002, component 1)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00238224 (tol = 0.002, component 1)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0126724 (tol = 0.002, component 1)print(steppedlarvareatFE)## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 12 -6182.0 12388
## (1 | TrtTankID) 1 11 -6182.0 12386 0.000 1 0.9977038
## (1 | AccTankID) 2 10 -6182.5 12385 1.085 1 0.2975686
## (1 | CrossID) 3 9 -6183.8 12386 2.571 1 0.1088464
## (1 | FemaleID) 0 8 -6202.1 12420 36.662 1 1.405e-09 ***
## (1 | MaleID) 0 8 -6190.8 12398 13.997 1 0.0001831 ***
## (1 | JarSeatable) 0 8 -6185.5 12387 3.418 1 0.0644882 .
## (1 | JarID) 0 8 -6294.2 12604 220.780 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## ParentTrt:JarTrt 0 26.811 26.811 1 240.66 5.3644 0.02139 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## sqrt(LarvaeAreaum2) ~ ParentTrt + JarTrt + (1 | FemaleID) + (1 |
## MaleID) + (1 | JarSeatable) + (1 | JarID) + ParentTrt:JarTrt#final model chosen by the lme4 step function: y~ParentTrt*JarTrt+(1|FemaleID)+(1|MaleID)+(1|JarID)+(1|JarSeatable)
larvareafintFE<- lmer(sqrt(LarvaeAreaum2)~ParentTrt*JarTrt+(1|FemaleID)+(1|MaleID)+(1|JarID)+(1|JarSeatable), data=LarvaeDat)
plot(larvareafintFE) #seems to meet assumption of linearity and homoscedasticity
qqnorm(resid(larvareafintFE))
qqline(resid(larvareafintFE)) #this seems to have improved it
acf(resid(larvareafintFE)) #this looks fine to me. 
summary(larvareafintFE)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: sqrt(LarvaeAreaum2) ~ ParentTrt * JarTrt + (1 | FemaleID) + (1 |
## MaleID) + (1 | JarID) + (1 | JarSeatable)
## Data: LarvaeDat
##
## REML criterion at convergence: 12367.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.4562 -0.5434 0.0596 0.6310 5.1810
##
## Random effects:
## Groups Name Variance Std.Dev.
## JarID (Intercept) 1.17016 1.0817
## FemaleID (Intercept) 0.55120 0.7424
## MaleID (Intercept) 0.26425 0.5141
## JarSeatable (Intercept) 0.05998 0.2449
## Residual 4.99788 2.2356
## Number of obs: 2698, groups:
## JarID, 270; FemaleID, 15; MaleID, 11; JarSeatable, 3
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 62.2179 0.4069 18.4716 152.903 <2e-16 ***
## ParentTrt2800 0.1959 0.5433 19.2917 0.361 0.7223
## JarTrtExposed -9.5031 0.2156 240.6370 -44.077 <2e-16 ***
## ParentTrt2800:JarTrtExposed 0.7307 0.3155 240.6577 2.316 0.0214 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) PrT2800 JrTrtE
## PrntTrt2800 -0.659
## JarTrtExpsd -0.269 0.201
## PrT2800:JTE 0.184 -0.292 -0.683Larvae Area Visualization
#plot the data using the means for each jar
ggplot(LarvaeDat,
aes(fill = WCa_pHDIC, y = LarvaeAreaum2, x = JarOmegaCalcite)) +
geom_point(aes(shape= Block, colour=WCa_pHDIC)) +
ggtitle("Larvae area by larvae treatment")+
ylab("Larvae Area (um2)")+
theme_classic()
#get means and SE of larvae areas
meanarea<- aggregate(LarvaeAreaum2~ParentTrt*JarTrt, data=LarvaeDat, FUN=mean)
searea<- aggregate(LarvaeAreaum2~ParentTrt*JarTrt, data=LarvaeDat, FUN=se)
#use a bargraph to plot the mean Larvae area +/- SE
ggplot(data = meanarea, aes(x = ParentTrt, y = LarvaeAreaum2, fill=JarTrt)) +
geom_bar(stat="identity", aes(fill=JarTrt), position="dodge")+
labs(x = "Parent Treatment", y = "Larvae Area (um2)") +
scale_fill_manual(name = "Larvae Treatment", labels = c("Control", "OA"), values = c("skyblue2", "red2")) +
geom_errorbar(aes(ymin= meanarea$LarvaeAreaum2-searea$LarvaeAreaum2, ymax=meanarea$LarvaeAreaum2+searea$LarvaeAreaum2),width=0.2, position=position_dodge(.9))+
theme_classic()## Warning: Use of `meanarea$LarvaeAreaum2` is discouraged. Use `LarvaeAreaum2`
## instead.
## Warning: Use of `meanarea$LarvaeAreaum2` is discouraged. Use `LarvaeAreaum2`
## instead.
#now plot the difference between parental treatments for larvae area
areadiff <- data.frame(tapply(LarvaeDat$LarvaeAreaum2, list(LarvaeDat$CrossID, LarvaeDat$JarTrt), mean)) %>%
rownames_to_column("CrossID") %>%
mutate(CrossTemp = CrossID) %>%
mutate(ParentTrt = ifelse(startsWith(CrossID, "E"), "Exposed", "Control")) %>%
separate(CrossTemp, c("FemaleID", "MaleID", "BlockID", "Fourth"), sep="_")%>%
mutate(Diff = Exposed - Control)
areadiff$CrossID<- as.factor(areadiff$CrossID)
areadiff$FemaleID<-as.factor(areadiff$FemaleID)
areadiff$MaleID<-as.factor(areadiff$MaleID)
#look at cross differences by Male and Female ID
ggplot(data = areadiff, aes(x = MaleID, y = Diff)) +
geom_point(aes(colour=FemaleID)) +
labs(x = "Male ID", y = "Change in mean shell area (um2) per family\nfrom control to OA conditions") +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), panel.background = element_blank(),
axis.line = element_line(colour = "black"), panel.border=element_rect(colour="black", fill=NA))+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5))
#get the means of the differences to visualize overall differences
meanareadiff<- aggregate(Diff~ParentTrt, data=areadiff, FUN=mean)
seareadiff<- aggregate(Diff~ParentTrt, data=areadiff, FUN=se)
ggplot(data = meanareadiff, aes(x = ParentTrt, y = Diff)) +
geom_bar(stat="identity", aes(fill="gray45"), position="dodge") +
labs(x = "Parental Treatment", y = "Change in mean shell area (um2) per family\nfrom control to OA conditions") +
scale_x_discrete(breaks = c("Control", "Exposed"), labels = c("Control", "OA")) +
scale_fill_manual(values = "gray45") +
geom_errorbar(aes(ymin= meanareadiff$Diff-seareadiff$Diff, ymax=meanareadiff$Diff+seareadiff$Diff),width=0.2, position=position_dodge(.9))+
theme(legend.position = "none", panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), panel.background = element_blank(),
axis.line = element_line(colour = "black"))## Warning: Use of `meanareadiff$Diff` is discouraged. Use `Diff` instead.## Warning: Use of `meanareadiff$Diff` is discouraged. Use `Diff` instead.
diffareattest<- t.test(Diff~ParentTrt, data=areadiff)
diffareattest #mean shell area difference is not significantly higher for OA parent treatment##
## Welch Two Sample t-test
##
## data: Diff by ParentTrt
## t = -1.7272, df = 39.691, p-value = 0.09192
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -168.53796 13.23648
## sample estimates:
## mean in group Control mean in group Exposed
## -1093.243 -1015.592Characterize the larvae by the ratio of major and minor axis
boxplot(MajMinRat~JarTrt*ParentTrt, data=LarvaeDat)
larvrat<-lmer(MajMinRat~WCa_pHDIC*JarOmegaCalcite+ (1|FemaleID)+(1|MaleID)+(1|JarSeatable)+(1|JarID)+(1|AccTankID) +(1|CrossID), data=LarvaeDat)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00219325 (tol = 0.002, component 1)larvratstep<- step(larvrat)## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingularprint(larvratstep)## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 11 4750.1 -9478.2
## (1 | AccTankID) 1 10 4750.1 -9480.1 0.047 1 0.8285
## (1 | JarSeatable) 2 9 4750.0 -9482.0 0.119 1 0.7298
## (1 | MaleID) 3 8 4749.8 -9483.6 0.395 1 0.5299
## (1 | CrossID) 4 7 4749.0 -9483.9 1.693 1 0.1932
## (1 | FemaleID) 0 6 4738.3 -9464.6 21.363 1 3.801e-06 ***
## (1 | JarID) 0 6 4731.8 -9451.6 34.298 1 4.728e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value
## WCa_pHDIC:JarOmegaCalcite 0 0.0068918 0.0068918 1 252.23 4.2751
## Pr(>F)
## WCa_pHDIC:JarOmegaCalcite 0.03969 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## MajMinRat ~ WCa_pHDIC + JarOmegaCalcite + (1 | FemaleID) + (1 |
## JarID) + WCa_pHDIC:JarOmegaCalcitelarvratfin<- lmer(MajMinRat~WCa_pHDIC*JarOmegaCalcite +(1|FemaleID)+(1|JarID), data=LarvaeDat)
#check the model to see that it meets assumptions
plot(larvratfin) #seems to meet assumption of linearity and homoscedasticity
qqnorm(resid(larvratfin))
qqline(resid(larvratfin)) #has relatively long tails, doesn't seem to meet assumption of normality. try a transformation
acf(resid(larvratfin)) #this looks fine.
summary(larvratfin) #significant effect of parent trt and significant interaction between parent treatment and larval treatment## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: MajMinRat ~ WCa_pHDIC * JarOmegaCalcite + (1 | FemaleID) + (1 |
## JarID)
## Data: LarvaeDat
##
## REML criterion at convergence: -9497.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.6168 -0.5759 -0.0263 0.5463 4.7286
##
## Random effects:
## Groups Name Variance Std.Dev.
## JarID (Intercept) 1.078e-04 0.010382
## FemaleID (Intercept) 5.034e-05 0.007095
## Residual 1.612e-03 0.040151
## Number of obs: 2698, groups: JarID, 270; FemaleID, 15
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.171609 0.005785 31.004778 202.522 <2e-16 ***
## WCa_pHDIC -0.013115 0.006313 30.846657 -2.077 0.0462 *
## JarOmegaCalcite -0.008297 0.005224 252.351564 -1.588 0.1135
## WCa_pHDIC:JarOmegaCalcite 0.011805 0.005710 252.231250 2.068 0.0397 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) WC_pHDIC JrOmgC
## WCa_pHDIC -0.893
## JarOmegClct -0.598 0.534
## WC_HDIC:JOC 0.533 -0.597 -0.893#try transforming the data
larvratt<-lmer(MajMinRat^1/3~WCa_pHDIC*JarOmegaCalcite+ (1|FemaleID)+(1|MaleID)+(1|JarSeatable)+(1|JarID)+(1|AccTankID) +(1|CrossID), data=LarvaeDat)
larvratstept<- step(larvratt)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00263927 (tol = 0.002, component 1)## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingularprint(larvratstept)## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 11 7709.8 -15398
## (1 | AccTankID) 1 10 7709.7 -15400 0.047 1 0.8285
## (1 | JarSeatable) 2 9 7709.7 -15401 0.119 1 0.7298
## (1 | MaleID) 3 8 7709.5 -15403 0.395 1 0.5299
## (1 | CrossID) 4 7 7708.6 -15403 1.693 1 0.1932
## (1 | FemaleID) 0 6 7697.9 -15384 21.363 1 3.801e-06 ***
## (1 | JarID) 0 6 7691.5 -15371 34.298 1 4.728e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value
## WCa_pHDIC:JarOmegaCalcite 0 0.00076575 0.00076575 1 252.23 4.2751
## Pr(>F)
## WCa_pHDIC:JarOmegaCalcite 0.03969 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## MajMinRat^1/3 ~ WCa_pHDIC + JarOmegaCalcite + (1 | FemaleID) +
## (1 | JarID) + WCa_pHDIC:JarOmegaCalcitelarvratfint<- lmer((MajMinRat^1/3)~WCa_pHDIC*JarOmegaCalcite +(1|FemaleID)+(1|JarID), data=LarvaeDat)
#check the model to see that it meets assumptions
plot(larvratfint) #seems to meet assumption of linearity and homoscedasticity
qqnorm(resid(larvratfint))
qqline(resid(larvratfint)) #has relatively long tails, doesn't seem to meet assumption of normality despite transformation
summary(larvratfint) #significant effect of parent trt and significant interaction between parent treatment and larval treatment## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: (MajMinRat^1/3) ~ WCa_pHDIC * JarOmegaCalcite + (1 | FemaleID) +
## (1 | JarID)
## Data: LarvaeDat
##
## REML criterion at convergence: -15417.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.6168 -0.5759 -0.0263 0.5463 4.7286
##
## Random effects:
## Groups Name Variance Std.Dev.
## JarID (Intercept) 1.198e-05 0.003461
## FemaleID (Intercept) 5.593e-06 0.002365
## Residual 1.791e-04 0.013384
## Number of obs: 2698, groups: JarID, 270; FemaleID, 15
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.390536 0.001928 31.004777 202.522 <2e-16 ***
## WCa_pHDIC -0.004372 0.002104 30.846656 -2.077 0.0462 *
## JarOmegaCalcite -0.002766 0.001741 252.351548 -1.588 0.1135
## WCa_pHDIC:JarOmegaCalcite 0.003935 0.001903 252.231235 2.068 0.0397 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) WC_pHDIC JrOmgC
## WCa_pHDIC -0.893
## JarOmegClct -0.598 0.534
## WC_HDIC:JOC 0.533 -0.597 -0.893#I have tried transforming with reciprocal, square root, square, and cube root none of them have helped with normality. No matter what though, still get a significant effect of parent treatment and significant interactionAnalyze the larvae by the ratio of major and minor axis with parent as a covariate and larvae as a fixed effect
larvratJF<-lmer(MajMinRat~WCa_pHDIC*JarTrt+ (1|FemaleID)+(1|MaleID)+(1|JarSeatable)+(1|JarID)+(1|AccTankID) +(1|CrossID), data=LarvaeDat)
larvratstepJF<- step(larvratJF)## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingularprint(larvratstepJF)## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 11 4749.8 -9477.5
## (1 | AccTankID) 1 10 4749.7 -9479.5 0.046 1 0.8293
## (1 | JarSeatable) 2 9 4749.7 -9481.4 0.114 1 0.7355
## (1 | MaleID) 3 8 4749.5 -9483.0 0.397 1 0.5289
## (1 | CrossID) 4 7 4748.6 -9483.3 1.705 1 0.1917
## (1 | FemaleID) 0 6 4737.9 -9463.9 21.416 1 3.697e-06 ***
## (1 | JarID) 0 6 4731.5 -9451.0 34.285 1 4.760e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## WCa_pHDIC:JarTrt 0 0.0069002 0.0069002 1 252.22 4.2804 0.03957 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## MajMinRat ~ WCa_pHDIC + JarTrt + (1 | FemaleID) + (1 | JarID) +
## WCa_pHDIC:JarTrtlarvratfinJF<- lmer(MajMinRat~WCa_pHDIC*JarTrt +(1|FemaleID)+(1|JarID), data=LarvaeDat)
#check the model to see that it meets assumptions
plot(larvratfinJF) #seems to meet assumption of linearity and homoscedasticity
qqnorm(resid(larvratfinJF))
qqline(resid(larvratfinJF)) #has relatively long tails, doesn't seem to meet assumption of normality. try a transformation
acf(resid(larvratfinJF)) #this looks fine.
summary(larvratfinJF)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: MajMinRat ~ WCa_pHDIC * JarTrt + (1 | FemaleID) + (1 | JarID)
## Data: LarvaeDat
##
## REML criterion at convergence: -9497.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.6172 -0.5756 -0.0274 0.5463 4.7283
##
## Random effects:
## Groups Name Variance Std.Dev.
## JarID (Intercept) 1.078e-04 0.010381
## FemaleID (Intercept) 5.043e-05 0.007101
## Residual 1.612e-03 0.040151
## Number of obs: 2698, groups: JarID, 270; FemaleID, 15
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.163e+00 5.144e-03 1.958e+01 226.008 <2e-16 ***
## WCa_pHDIC -2.818e-04 5.627e-03 1.966e+01 -0.050 0.9606
## JarTrtExposed 7.059e-03 4.439e-03 2.523e+02 1.590 0.1130
## WCa_pHDIC:JarTrtExposed -1.003e-02 4.850e-03 2.522e+02 -2.069 0.0396 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) WC_pHDIC JrTrtE
## WCa_pHDIC -0.893
## JarTrtExpsd -0.432 0.386
## WC_HDIC:JTE 0.386 -0.434 -0.893Analyze the larvae by the ratio of major and minor axis with parent and larvae as fixed effects
larvratFE<-lmer(MajMinRat~ParentTrt*JarTrt+ (1|FemaleID)+(1|MaleID)+(1|JarSeatable)+(1|JarID)+(1|TrtTankID)+(1|AccTankID) +(1|CrossID), data=LarvaeDat)
larvratstepFE<- step(larvratFE)## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular## Warning: Model failed to converge with 1 negative eigenvalue: -2.9e+01## boundary (singular) fit: see ?isSingularprint(larvratstepFE)## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 12 4749.5 -9475.1
## (1 | TrtTankID) 1 11 4749.5 -9477.1 0.000 1 0.9991
## (1 | AccTankID) 2 10 4749.5 -9479.0 0.047 1 0.8282
## (1 | JarSeatable) 3 9 4749.5 -9480.9 0.116 1 0.7331
## (1 | MaleID) 4 8 4749.3 -9482.5 0.392 1 0.5312
## (1 | CrossID) 5 7 4748.4 -9482.8 1.730 1 0.1884
## (1 | FemaleID) 0 6 4737.4 -9462.8 22.000 1 2.727e-06 ***
## (1 | JarID) 0 6 4731.4 -9450.8 34.009 1 5.486e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## ParentTrt:JarTrt 0 0.0076797 0.0076797 1 252.22 4.7639 0.02999 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## MajMinRat ~ ParentTrt + JarTrt + (1 | FemaleID) + (1 | JarID) +
## ParentTrt:JarTrtlarvratfinFE<- lmer(MajMinRat~ParentTrt*JarTrt +(1|FemaleID)+(1|JarID), data=LarvaeDat)
#check the model to see that it meets assumptions
plot(larvratfinFE) #seems to meet assumption of linearity and homoscedasticity
qqnorm(resid(larvratfinFE))
qqline(resid(larvratfinFE)) #has relatively long tails, doesn't seem to meet assumption of normality
acf(resid(larvratfinFE)) #this looks fine.
summary(larvratfinFE)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: MajMinRat ~ ParentTrt * JarTrt + (1 | FemaleID) + (1 | JarID)
## Data: LarvaeDat
##
## REML criterion at convergence: -9496.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.6210 -0.5758 -0.0245 0.5460 4.7281
##
## Random effects:
## Groups Name Variance Std.Dev.
## JarID (Intercept) 0.0001072 0.010356
## FemaleID (Intercept) 0.0000514 0.007169
## Residual 0.0016121 0.040151
## Number of obs: 2698, groups: JarID, 270; FemaleID, 15
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.163e+00 3.195e-03 1.966e+01 363.863 <2e-16
## ParentTrt2800 -3.787e-04 4.670e-03 1.955e+01 -0.081 0.9362
## JarTrtExposed -5.217e-03 2.733e-03 2.521e+02 -1.909 0.0574
## ParentTrt2800:JarTrtExposed 8.730e-03 4.000e-03 2.522e+02 2.183 0.0300
##
## (Intercept) ***
## ParentTrt2800
## JarTrtExposed .
## ParentTrt2800:JarTrtExposed *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) PrT2800 JrTrtE
## PrntTrt2800 -0.684
## JarTrtExpsd -0.434 0.297
## PrT2800:JTE 0.296 -0.431 -0.683Visualize the ratio of the axes
#get the mean ratio for the treatments to plot
meanrat<- ddply(LarvaeDat, .(JarTrt, ParentTrt), numcolwise(mean, na.rm=T))
serat<- ddply(LarvaeDat, .(JarTrt, ParentTrt), numcolwise(se, na.rm=T))
ggplot(LarvaeDat, aes(x=MajMinRat, color=JarTrt))+
geom_histogram(alpha=0.5, position="identity") +
facet_grid(~ParentTrt)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(data = meanrat, aes(x = ParentTrt, y = MajMinRat, fill=JarTrt)) +
geom_bar(stat="identity", aes(fill=JarTrt), position="dodge")+
labs(x = "Parent Treatment", y = "Ratio of Major to Minor Axis") +
scale_fill_manual(name = "Larvae Treatment", labels = c("Control", "OA Exposed"), values = c("skyblue2", "red2")) +
geom_errorbar(aes(ymin= meanrat$MajMinRat-serat$MajMinRat, ymax=meanrat$MajMinRat+serat$MajMinRat),width=0.2, position=position_dodge(.9))+
theme_classic()## Warning: Use of `meanrat$MajMinRat` is discouraged. Use `MajMinRat` instead.
## Warning: Use of `meanrat$MajMinRat` is discouraged. Use `MajMinRat` instead.
Look at the ratio of perimeter to diameter. For a circle, this would be equal to pi
PerDiRat<- lmer(PerimDiamRat~WCa_pHDIC*JarOmegaCalcite+(1|FemaleID)+(1|MaleID)+(1|JarSeatable)+(1|JarID)+(1|AccTankID) +(1|CrossID), data=LarvaeDat)## boundary (singular) fit: see ?isSingularPerDiStep<- step(PerDiRat)## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular## Warning: Model failed to converge with 1 negative eigenvalue: -1.2e+01## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingularprint(PerDiStep)## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 11 4346.6 -8671.2
## (1 | MaleID) 1 10 4346.6 -8673.2 0.000 1 1
## (1 | AccTankID) 2 9 4346.6 -8675.2 0.000 1 1
## (1 | JarSeatable) 3 8 4346.6 -8677.2 0.000 1 1
## (1 | CrossID) 4 7 4346.6 -8679.2 0.000 1 1
## (1 | FemaleID) 0 6 4337.6 -8663.2 17.974 1 2.239e-05 ***
## (1 | JarID) 0 6 4333.8 -8655.5 25.700 1 3.989e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value
## WCa_pHDIC:JarOmegaCalcite 0 0.011272 0.011272 1 253 5.1505
## Pr(>F)
## WCa_pHDIC:JarOmegaCalcite 0.02408 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## PerimDiamRat ~ WCa_pHDIC + JarOmegaCalcite + (1 | FemaleID) +
## (1 | JarID) + WCa_pHDIC:JarOmegaCalcite#final model chosen PerimDiamRat ~ WCa_pHDIC * JarOmegaCalcite + (1 | FemaleID) + (1 | JarID)
PerDiRatFin<- lmer(PerimDiamRat ~ WCa_pHDIC * JarOmegaCalcite + (1 | FemaleID) + (1 | JarID), data=LarvaeDat)
#check the model to see that it meets assumptions
plot(PerDiRatFin) #seems to meet assumption of linearity and homoscedasticity
qqnorm(resid(PerDiRatFin))
qqline(resid(PerDiRatFin)) #this looks fine
acf(resid(PerDiRatFin))
summary(PerDiRatFin) #significant effect of parent treatment, jar, and interaction, not an additive effect though. ## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: PerimDiamRat ~ WCa_pHDIC * JarOmegaCalcite + (1 | FemaleID) +
## (1 | JarID)
## Data: LarvaeDat
##
## REML criterion at convergence: -8693.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.1373 -0.5713 0.0461 0.6356 3.1915
##
## Random effects:
## Groups Name Variance Std.Dev.
## JarID (Intercept) 1.231e-04 0.011094
## FemaleID (Intercept) 5.609e-05 0.007489
## Residual 2.188e-03 0.046781
## Number of obs: 2698, groups: JarID, 270; FemaleID, 15
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.084977 0.006320 33.245328 488.133 < 2e-16 ***
## WCa_pHDIC 0.018379 0.006896 33.063404 2.665 0.0118 *
## JarOmegaCalcite 0.049719 0.005889 253.121060 8.442 2.44e-15 ***
## WCa_pHDIC:JarOmegaCalcite -0.014609 0.006437 252.997219 -2.269 0.0241 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) WC_pHDIC JrOmgC
## WCa_pHDIC -0.893
## JarOmegClct -0.617 0.551
## WC_HDIC:JOC 0.550 -0.616 -0.893#larvae in control treatments from both parents are more circular. Larvae are shorter and squatter from OA treatments. If they were elongated, the value would be greater than pi. Look at the ratio of perimeter to diameter with Parent as a covariate and JarTrt as a fixed effect. For a circle, this would be equal to pi
PerDiRatJF<- lmer(PerimDiamRat~WCa_pHDIC*JarTrt+(1|FemaleID)+(1|MaleID)+(1|JarSeatable)+(1|JarID)+(1|AccTankID) +(1|CrossID), data=LarvaeDat)## boundary (singular) fit: see ?isSingularPerDiStepJF<- step(PerDiRatJF)## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingularprint(PerDiStepJF)## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 11 4346.7 -8671.4
## (1 | AccTankID) 1 10 4346.7 -8673.4 0.000 1 1
## (1 | MaleID) 2 9 4346.7 -8675.4 0.000 1 1
## (1 | JarSeatable) 3 8 4346.7 -8677.4 0.000 1 1
## (1 | CrossID) 4 7 4346.7 -8679.4 0.000 1 1
## (1 | FemaleID) 0 6 4337.6 -8663.3 18.124 1 2.07e-05 ***
## (1 | JarID) 0 6 4334.0 -8656.1 25.295 1 4.92e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## WCa_pHDIC:JarTrt 0 0.011719 0.011719 1 252.98 5.3549 0.02147 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## PerimDiamRat ~ WCa_pHDIC + JarTrt + (1 | FemaleID) + (1 | JarID) +
## WCa_pHDIC:JarTrt#final model chosen PerimDiamRat ~ WCa_pHDIC * JarTrt + (1 | FemaleID) + (1 | JarID)
PerDiRatFinJF<- lmer(PerimDiamRat ~ WCa_pHDIC * JarTrt + (1 | FemaleID) + (1 | JarID), data=LarvaeDat)
#check the model to see that it meets assumptions
plot(PerDiRatFinJF) #seems to meet assumption of linearity and homoscedasticity
qqnorm(resid(PerDiRatFinJF))
qqline(resid(PerDiRatFinJF)) #this looks fine
acf(resid(PerDiRatFinJF))
summary(PerDiRatFinJF)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: PerimDiamRat ~ WCa_pHDIC * JarTrt + (1 | FemaleID) + (1 | JarID)
## Data: LarvaeDat
##
## REML criterion at convergence: -8693.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.1399 -0.5686 0.0462 0.6360 3.1972
##
## Random effects:
## Groups Name Variance Std.Dev.
## JarID (Intercept) 1.219e-04 0.01104
## FemaleID (Intercept) 5.625e-05 0.00750
## Residual 2.188e-03 0.04678
## Number of obs: 2698, groups: JarID, 270; FemaleID, 15
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.139138 0.005569 20.307544 563.724 < 2e-16 ***
## WCa_pHDIC 0.002370 0.006091 20.396908 0.389 0.7012
## JarTrtExposed -0.042459 0.004997 253.096458 -8.498 1.68e-15 ***
## WCa_pHDIC:JarTrtExposed 0.012634 0.005460 252.978240 2.314 0.0215 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) WC_pHDIC JrTrtE
## WCa_pHDIC -0.893
## JarTrtExpsd -0.449 0.401
## WC_HDIC:JTE 0.402 -0.451 -0.893Look at the ratio of perimeter to diameter with Parent and JarTrt as a fixed effects
PerDiRatFE<- lmer(PerimDiamRat~ParentTrt*JarTrt+(1|FemaleID)+(1|MaleID)+(1|JarSeatable)+(1|JarID)+(1|AccTankID) +(1|CrossID)+(1|TrtTankID), data=LarvaeDat)## boundary (singular) fit: see ?isSingularPerDiStepFE<- step(PerDiRatFE)## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingularprint(PerDiStepFE)## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 12 4346.4 -8668.8
## (1 | MaleID) 1 11 4346.4 -8670.8 0.000 1 1
## (1 | AccTankID) 2 10 4346.4 -8672.8 0.000 1 1
## (1 | TrtTankID) 3 9 4346.4 -8674.8 0.000 1 1
## (1 | CrossID) 4 8 4346.4 -8676.8 0.000 1 1
## (1 | JarSeatable) 5 7 4346.4 -8678.8 0.000 1 1
## (1 | FemaleID) 0 6 4336.9 -8661.7 19.087 1 1.249e-05 ***
## (1 | JarID) 0 6 4333.9 -8655.8 25.024 1 5.663e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## ParentTrt:JarTrt 0 0.012977 0.012977 1 252.99 5.9295 0.01558 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## PerimDiamRat ~ ParentTrt + JarTrt + (1 | FemaleID) + (1 | JarID) +
## ParentTrt:JarTrt#final model chosen PerimDiamRat ~ WCa_pHDIC * JarTrt + (1 | FemaleID) + (1 | JarID)
PerDiRatFinFE<- lmer(PerimDiamRat ~ ParentTrt * JarTrt + (1 | FemaleID) + (1 | JarID), data=LarvaeDat)
#check the model to see that it meets assumptions
plot(PerDiRatFinFE) #seems to meet assumption of linearity and homoscedasticity
qqnorm(resid(PerDiRatFinFE))
qqline(resid(PerDiRatFinFE)) #this looks fine
acf(resid(PerDiRatFinFE))
summary(PerDiRatFinFE)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: PerimDiamRat ~ ParentTrt * JarTrt + (1 | FemaleID) + (1 | JarID)
## Data: LarvaeDat
##
## REML criterion at convergence: -8692.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.1404 -0.5705 0.0466 0.6357 3.2012
##
## Random effects:
## Groups Name Variance Std.Dev.
## JarID (Intercept) 1.212e-04 0.01101
## FemaleID (Intercept) 5.837e-05 0.00764
## Residual 2.188e-03 0.04678
## Number of obs: 2698, groups: JarID, 270; FemaleID, 15
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 3.141620 0.003477 20.281561 903.594 < 2e-16
## ParentTrt2800 -0.001170 0.005082 20.162026 -0.230 0.8202
## JarTrtExposed -0.027015 0.003075 252.818361 -8.785 2.43e-16
## ParentTrt2800:JarTrtExposed -0.010961 0.004501 252.989564 -2.435 0.0156
##
## (Intercept) ***
## ParentTrt2800
## JarTrtExposed ***
## ParentTrt2800:JarTrtExposed *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) PrT2800 JrTrtE
## PrntTrt2800 -0.684
## JarTrtExpsd -0.448 0.307
## PrT2800:JTE 0.306 -0.446 -0.683Visualize Perimeter to Diameter ratio
ggplot(LarvaeDat, aes(x=PerimDiamRat, color=JarTrt))+
geom_histogram(alpha=0.5, position="identity") +
facet_grid(~ParentTrt)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
meanperimdiam<- aggregate(PerimDiamRat~ParentTrt*JarTrt, data=LarvaeDat, FUN=mean)
seperimdiam<- aggregate(PerimDiamRat~ParentTrt*JarTrt, data=LarvaeDat, FUN=se)
ggplot(data = meanperimdiam, aes(x = ParentTrt, y = PerimDiamRat, fill=JarTrt)) +
geom_bar(stat="identity", aes(fill=JarTrt), position="dodge")+
labs(x = "Parent Treatment", y = "Perimeter to Diameter Ratio") +
scale_fill_manual(name = "Larvae Treatment", labels = c("Control", "OA"), values = c("skyblue2", "red2")) +
geom_errorbar(aes(ymin= PerimDiamRat-seperimdiam$PerimDiamRat, ymax=PerimDiamRat+seperimdiam$PerimDiamRat),width=0.2, position=position_dodge(.9))+
geom_abline(slope=0, intercept=pi,colour="orange")+
theme_classic()
Look at eccentricity
#look at the eccentricity data
hist(LarvaeDat$LarvaeEccentricity)#data apprear to be left skewed. 
#get the mean eccentricity for the treatments to plot
meaneccen<- aggregate(LarvaeEccentricity~ParentTrt*JarTrt, data=LarvaeDat, FUN=mean)
seeccen<- aggregate(LarvaeEccentricity~ParentTrt*JarTrt, data=LarvaeDat, FUN=se)
ggplot(data = meaneccen, aes(x = ParentTrt, y = LarvaeEccentricity, fill=JarTrt)) +
geom_bar(stat="identity", aes(fill=JarTrt), position="dodge")+
labs(x = "Parent Treatment", y = "Larvae Eccentricity") +
scale_fill_manual(name = "Larvae Treatment", labels = c("Control", "OA"), values = c("skyblue2", "red2")) +
geom_errorbar(aes(ymin= LarvaeEccentricity-seeccen$LarvaeEccentricity, ymax=LarvaeEccentricity+seeccen$LarvaeEccentricity),width=0.2, position=position_dodge(.9))+
theme_classic()
jarmean<-aggregate(LarvaeEccentricity~JarTrt, data=LarvaeDat, FUN=mean)
sejar<- aggregate(LarvaeEccentricity~JarTrt, data=LarvaeDat, FUN=se)
#plot just by JarTrt
ggplot(data = jarmean, aes(x = JarTrt, y = LarvaeEccentricity)) +
geom_bar(stat="identity", aes(fill="gray45"), position="dodge")+
labs(x = "Jar Treatment", y = "Larvae Eccentricity") +
scale_fill_manual(values = c("gray45")) +
geom_errorbar(aes(ymin= LarvaeEccentricity-sejar$LarvaeEccentricity, ymax=LarvaeEccentricity+sejar$LarvaeEccentricity),width=0.2, position=position_dodge(.9))+
theme_classic()
larveccen<-lmer(LarvaeEccentricity~WCa_pHDIC*JarOmegaCalcite+ (1|FemaleID)+(1|MaleID)+(1|CrossID)+(1|JarSeatable)+(1|JarID)+(1|AccTankID), data=LarvaeDat)
steppedlarveccen<- step(larveccen)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00337866 (tol = 0.002, component 1)## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingularprint(steppedlarveccen)## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 11 3928.2 -7834.3
## (1 | AccTankID) 1 10 3928.2 -7836.3 0.000 1 1.0000
## (1 | MaleID) 2 9 3928.2 -7838.3 0.003 1 0.9555
## (1 | JarSeatable) 3 8 3928.1 -7840.2 0.083 1 0.7732
## (1 | CrossID) 4 7 3927.4 -7840.8 1.446 1 0.2292
## (1 | FemaleID) 0 6 3917.7 -7823.5 19.345 1 1.091e-05 ***
## (1 | JarID) 0 6 3906.9 -7801.8 40.998 1 1.524e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value
## WCa_pHDIC:JarOmegaCalcite 1 0.010813 0.010813 1 252.181 3.6597
## WCa_pHDIC 2 0.003051 0.003051 1 13.003 1.0327
## JarOmegaCalcite 0 0.021731 0.021731 1 253.313 7.3556
## Pr(>F)
## WCa_pHDIC:JarOmegaCalcite 0.056875 .
## WCa_pHDIC 0.328064
## JarOmegaCalcite 0.007144 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## LarvaeEccentricity ~ JarOmegaCalcite + (1 | FemaleID) + (1 |
## JarID)#final model only had JarTrt, female ID and jarID
larveccenfin<- lmer(LarvaeEccentricity~JarOmegaCalcite+ (1|FemaleID)+(1|JarID), data=LarvaeDat)
#check the model to see that it meets assumptions
plot(larveccenfin) #seems to meet assumption of linearity and homoscedasticity
qqnorm(resid(larveccenfin))
qqline(resid(larveccenfin)) #has relatively long tails, doesn't meet assumption of normality. 
acf(resid(larveccenfin))
summary(larveccenfin) #significant effect of JarTrt## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: LarvaeEccentricity ~ JarOmegaCalcite + (1 | FemaleID) + (1 |
## JarID)
## Data: LarvaeDat
##
## REML criterion at convergence: -7866.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.8185 -0.4674 0.0592 0.5861 3.4958
##
## Random effects:
## Groups Name Variance Std.Dev.
## JarID (Intercept) 2.256e-04 0.015020
## FemaleID (Intercept) 8.923e-05 0.009446
## Residual 2.954e-03 0.054354
## Number of obs: 2698, groups: JarID, 270; FemaleID, 15
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 4.970e-01 3.540e-03 3.476e+01 140.388 < 2e-16 ***
## JarOmegaCalcite 8.879e-03 3.274e-03 2.533e+02 2.712 0.00714 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## JarOmegClct -0.609#try a transformation
larveccent<-lmer(LarvaeEccentricity^3~JarOmegaCalcite+ (1|FemaleID)+(1|JarID) , data=LarvaeDat)
#check the model to see that it meets assumptions
plot(larveccent) #seems to meet assumption of linearity and homoscedasticity
qqnorm(resid(larveccent))
qqline(resid(larveccent)) #has relatively long tails, doesn't meet assumption of normality. transformation helped, but it's not great
summary(larveccent) #no significant differences ## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: LarvaeEccentricity^3 ~ JarOmegaCalcite + (1 | FemaleID) + (1 |
## JarID)
## Data: LarvaeDat
##
## REML criterion at convergence: -9778.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2074 -0.5825 -0.0178 0.5619 4.4271
##
## Random effects:
## Groups Name Variance Std.Dev.
## JarID (Intercept) 1.006e-04 0.010028
## FemaleID (Intercept) 4.613e-05 0.006792
## Residual 1.459e-03 0.038201
## Number of obs: 2698, groups: JarID, 270; FemaleID, 15
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.311e-01 2.488e-03 3.317e+01 52.704 <2e-16 ***
## JarOmegaCalcite 1.134e-03 2.252e-03 2.534e+02 0.504 0.615
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## JarOmegClct -0.596Analyze eccentricity with parent as a covariate and jar as fixed effect
larveccenJF<-lmer(LarvaeEccentricity~WCa_pHDIC*JarTrt+ (1|FemaleID)+(1|MaleID)+(1|CrossID)+(1|JarSeatable)+(1|JarID)+(1|AccTankID), data=LarvaeDat)## boundary (singular) fit: see ?isSingular## Warning: Model failed to converge with 1 negative eigenvalue: -4.3e+01steppedlarveccenJF<- step(larveccenJF)## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular## Warning: Model failed to converge with 1 negative eigenvalue: -4.2e+01## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## unable to evaluate scaled gradient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge: degenerate Hessian with 1 negative eigenvalues## Warning: Model failed to converge with 1 negative eigenvalue: -2.5e+01## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## unable to evaluate scaled gradient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge: degenerate Hessian with 1 negative eigenvalues## Warning: Model failed to converge with 1 negative eigenvalue: -2.5e+01## boundary (singular) fit: see ?isSingularprint(steppedlarveccenJF)## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 11 3927.8 -7833.7
## (1 | AccTankID) 1 10 3927.8 -7835.7 0.000 1 1.0000
## (1 | MaleID) 2 9 3927.8 -7837.7 0.000 1 0.9995
## (1 | JarSeatable) 3 8 3927.8 -7839.6 0.087 1 0.7682
## (1 | CrossID) 4 7 3927.1 -7840.1 1.439 1 0.2304
## (1 | FemaleID) 0 6 3917.4 -7822.8 19.342 1 1.093e-05 ***
## (1 | JarID) 0 6 3906.6 -7801.1 41.002 1 1.521e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## WCa_pHDIC:JarTrt 1 0.0108336 0.0108336 1 252.166 3.6668 0.056637
## WCa_pHDIC 2 0.0030599 0.0030599 1 13.003 1.0357 0.327385
## JarTrt 0 0.0216891 0.0216891 1 253.301 7.3414 0.007198
##
## WCa_pHDIC:JarTrt .
## WCa_pHDIC
## JarTrt **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## LarvaeEccentricity ~ JarTrt + (1 | FemaleID) + (1 | JarID)#final model
larveccenfinJF<- lmer(LarvaeEccentricity~JarTrt+ (1|FemaleID)+(1|JarID), data=LarvaeDat)
#check the model to see that it meets assumptions
plot(larveccenfinJF) #seems to meet assumption of linearity and homoscedasticity
qqnorm(resid(larveccenfinJF))
qqline(resid(larveccenfinJF)) #has relatively long tails, doesn't meet assumption of normality. 
acf(resid(larveccenfinJF))
summary(larveccenfinJF)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: LarvaeEccentricity ~ JarTrt + (1 | FemaleID) + (1 | JarID)
## Data: LarvaeDat
##
## REML criterion at convergence: -7865.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.8185 -0.4705 0.0593 0.5863 3.4958
##
## Random effects:
## Groups Name Variance Std.Dev.
## JarID (Intercept) 2.257e-04 0.015023
## FemaleID (Intercept) 8.905e-05 0.009437
## Residual 2.954e-03 0.054354
## Number of obs: 2698, groups: JarID, 270; FemaleID, 15
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.506682 0.003135 21.768315 161.632 <2e-16 ***
## JarTrtExposed -0.007533 0.002780 253.301235 -2.709 0.0072 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## JarTrtExpsd -0.447Analyze eccentricity with parent and jar as fixed effects
larveccenFE<-lmer(LarvaeEccentricity~ParentTrt*JarTrt+ (1|FemaleID)+(1|MaleID)+(1|CrossID)+(1|JarSeatable)+(1|JarID)+(1|TrtTankID)+(1|AccTankID), data=LarvaeDat)## boundary (singular) fit: see ?isSingularsteppedlarveccenFE<- step(larveccenFE)## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular## Warning: Model failed to converge with 1 negative eigenvalue: -1.4e+00## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00459713 (tol = 0.002, component 1)## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingularprint(steppedlarveccenFE)## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 12 3927.5 -7831.1
## (1 | TrtTankID) 1 11 3927.5 -7833.1 0.000 1 1.0000
## (1 | AccTankID) 2 10 3927.5 -7835.1 0.003 1 0.9569
## (1 | MaleID) 3 9 3927.5 -7837.1 0.000 1 1.0000
## (1 | JarSeatable) 4 8 3927.5 -7839.0 0.088 1 0.7664
## (1 | CrossID) 5 7 3926.8 -7839.5 1.455 1 0.2277
## (1 | FemaleID) 0 6 3916.8 -7821.7 19.840 1 8.418e-06 ***
## (1 | JarID) 0 6 3906.4 -7800.8 40.785 1 1.699e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## ParentTrt:JarTrt 0 0.011839 0.011839 1 252.18 4.0073 0.04637 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## LarvaeEccentricity ~ ParentTrt + JarTrt + (1 | FemaleID) + (1 |
## JarID) + ParentTrt:JarTrt#final model
larveccenfinFE<- lmer(LarvaeEccentricity~ParentTrt*JarTrt+ (1|FemaleID)+(1|JarID), data=LarvaeDat)
#check the model to see that it meets assumptions
plot(larveccenfinFE) #seems to meet assumption of linearity and homoscedasticity
qqnorm(resid(larveccenfinFE))
qqline(resid(larveccenfinFE)) #has relatively long tails, doesn't meet assumption of normality. 
acf(resid(larveccenfinFE))
summary(larveccenfinFE)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: LarvaeEccentricity ~ ParentTrt * JarTrt + (1 | FemaleID) + (1 |
## JarID)
## Data: LarvaeDat
##
## REML criterion at convergence: -7853.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -6.7892 -0.4764 0.0642 0.5809 3.5249
##
## Random effects:
## Groups Name Variance Std.Dev.
## JarID (Intercept) 2.193e-04 0.014809
## FemaleID (Intercept) 9.099e-05 0.009539
## Residual 2.954e-03 0.054355
## Number of obs: 2698, groups: JarID, 270; FemaleID, 15
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 5.069e-01 4.316e-03 2.012e+01 117.439 < 2e-16
## ParentTrt2800 -3.682e-04 6.309e-03 2.000e+01 -0.058 0.954041
## JarTrtExposed -1.270e-02 3.784e-03 2.520e+02 -3.356 0.000913
## ParentTrt2800:JarTrtExposed 1.109e-02 5.538e-03 2.522e+02 2.002 0.046375
##
## (Intercept) ***
## ParentTrt2800
## JarTrtExposed ***
## ParentTrt2800:JarTrtExposed *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) PrT2800 JrTrtE
## PrntTrt2800 -0.684
## JarTrtExpsd -0.444 0.304
## PrT2800:JTE 0.304 -0.442 -0.683Test Larvae perimeter
ggplot(LarvaeDat, aes(x=LarvaePerimeterum, color=JarTrt))+
geom_histogram(alpha=0.5, position="identity") +
facet_grid(~ParentTrt)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
#test to see if there is a significant difference
larvperim<-lmer(LarvaePerimeterum~WCa_pHDIC*JarOmegaCalcite+ (1|FemaleID)+(1|MaleID)+(1|CrossID)+(1|JarSeatable)+(1|JarID)+(1|AccTankID), data=LarvaeDat)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0059559 (tol = 0.002, component 1)steppedlarvperim<- step(larvperim)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00511491 (tol = 0.002, component 1)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0027975 (tol = 0.002, component 1)print(steppedlarvperim)## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 11 -9769.7 19561
## (1 | AccTankID) 1 10 -9769.9 19560 0.459 1 0.4980403
## (1 | MaleID) 2 9 -9770.9 19560 1.951 1 0.1624814
## (1 | FemaleID) 0 8 -9773.5 19563 5.260 1 0.0218204 *
## (1 | CrossID) 0 8 -9777.3 19571 12.689 1 0.0003677 ***
## (1 | JarSeatable) 0 8 -9773.4 19563 4.965 1 0.0258594 *
## (1 | JarID) 0 8 -9869.4 19755 197.038 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value
## WCa_pHDIC:JarOmegaCalcite 0 573.37 573.37 1 221.4 8.0187
## Pr(>F)
## WCa_pHDIC:JarOmegaCalcite 0.005056 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## LarvaePerimeterum ~ WCa_pHDIC + JarOmegaCalcite + (1 | FemaleID) +
## (1 | CrossID) + (1 | JarSeatable) + (1 | JarID) + WCa_pHDIC:JarOmegaCalcite#final model chosen by the step function: LarvaePerimeterum~WCa_pHDIC*JarOmegaCalcite+ (1|FemaleID)+(1|CrossID)+(1|JarSeatable)+(1|JarID)
larvperimfin<- lmer(LarvaePerimeterum~WCa_pHDIC*JarOmegaCalcite+ (1|FemaleID)+(1|CrossID)+(1|JarSeatable)+(1|JarID), data=LarvaeDat)
#check the model to see that it meets assumptions
plot(larvperimfin) #seems to meet assumption of linearity and homoscedasticity
qqnorm(resid(larvperimfin))
qqline(resid(larvperimfin)) #has relatively long tails, doesn't meet assumption of normality. 
acf(resid(larvperimfin)) #this looks fine to me.
summary(larvperimfin) #should I try to transform? Main effects and interaction are significant. ## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: LarvaePerimeterum ~ WCa_pHDIC * JarOmegaCalcite + (1 | FemaleID) +
## (1 | CrossID) + (1 | JarSeatable) + (1 | JarID)
## Data: LarvaeDat
##
## REML criterion at convergence: 19541.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.4544 -0.5443 0.0484 0.6225 5.3183
##
## Random effects:
## Groups Name Variance Std.Dev.
## JarID (Intercept) 16.212 4.026
## CrossID (Intercept) 5.489 2.343
## FemaleID (Intercept) 5.599 2.366
## JarSeatable (Intercept) 1.107 1.052
## Residual 71.504 8.456
## Number of obs: 2698, groups:
## JarID, 270; CrossID, 45; FemaleID, 15; JarSeatable, 3
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 196.343 2.069 24.736 94.887 < 2e-16 ***
## WCa_pHDIC -5.511 2.160 23.829 -2.552 0.01755 *
## JarOmegaCalcite 37.384 1.539 221.250 24.284 < 2e-16 ***
## WCa_pHDIC:JarOmegaCalcite 4.767 1.683 221.395 2.832 0.00506 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) WC_pHDIC JrOmgC
## WCa_pHDIC -0.854
## JarOmegClct -0.493 0.460
## WC_HDIC:JOC 0.439 -0.514 -0.893Analyze larvae perimeter with adult as a covariate and jar as a fixed effect
larvperimJF<-lmer(LarvaePerimeterum~WCa_pHDIC*JarTrt+ (1|FemaleID)+(1|MaleID)+(1|CrossID)+(1|JarSeatable)+(1|JarID)+(1|AccTankID), data=LarvaeDat)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00324995 (tol = 0.002, component 1)steppedlarvperimJF<- step(larvperimJF)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00257919 (tol = 0.002, component 1)print(steppedlarvperimJF)## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 11 -9767.5 19557
## (1 | AccTankID) 1 10 -9767.7 19555 0.454 1 0.5003185
## (1 | CrossID) 2 9 -9769.0 19556 2.586 1 0.1077988
## (1 | FemaleID) 0 8 -9788.4 19593 38.693 1 4.959e-10 ***
## (1 | MaleID) 0 8 -9776.4 19569 14.715 1 0.0001251 ***
## (1 | JarSeatable) 0 8 -9772.4 19561 6.684 1 0.0097280 **
## (1 | JarID) 0 8 -9880.6 19777 223.234 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## WCa_pHDIC:JarTrt 0 505.5 505.5 1 240.59 7.0697 0.008366 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## LarvaePerimeterum ~ WCa_pHDIC + JarTrt + (1 | FemaleID) + (1 |
## MaleID) + (1 | JarSeatable) + (1 | JarID) + WCa_pHDIC:JarTrt#final model chosen by the step function: LarvaePerimeterum~WCa_pHDIC*JarTrt+ (1|FemaleID)+(1|MaleID)+(1|JarSeatable)+(1|JarID)
larvperimfinJF<- lmer(LarvaePerimeterum~WCa_pHDIC*JarTrt+ (1|FemaleID)+(1|MaleID)+(1|JarSeatable)+(1|JarID), data=LarvaeDat)
#check the model to see that it meets assumptions
plot(larvperimfinJF) #seems to meet assumption of linearity and homoscedasticity
qqnorm(resid(larvperimfinJF))
qqline(resid(larvperimfinJF)) #has relatively long tails, doesn't meet assumption of normality. 
acf(resid(larvperimfinJF)) #this looks fine to me.
summary(larvperimfinJF)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: LarvaePerimeterum ~ WCa_pHDIC * JarTrt + (1 | FemaleID) + (1 |
## MaleID) + (1 | JarSeatable) + (1 | JarID)
## Data: LarvaeDat
##
## REML criterion at convergence: 19538
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.3527 -0.5461 0.0479 0.6246 5.2579
##
## Random effects:
## Groups Name Variance Std.Dev.
## JarID (Intercept) 16.880 4.108
## FemaleID (Intercept) 8.381 2.895
## MaleID (Intercept) 3.937 1.984
## JarSeatable (Intercept) 1.432 1.197
## Residual 71.503 8.456
## Number of obs: 2698, groups:
## JarID, 270; FemaleID, 15; MaleID, 11; JarSeatable, 3
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 237.0684 2.3822 21.0475 99.516 < 2e-16 ***
## WCa_pHDIC -0.4725 2.5431 19.3191 -0.186 0.85453
## JarTrtExposed -31.8749 1.3265 240.6127 -24.029 < 2e-16 ***
## WCa_pHDIC:JarTrtExposed -3.8550 1.4499 240.5926 -2.659 0.00837 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) WC_pHDIC JrTrtE
## WCa_pHDIC -0.849
## JarTrtExpsd -0.279 0.255
## WC_HDIC:JTE 0.249 -0.287 -0.893Analyze larvae perimeter with adult and jar as a fixed effects
larvperimFE<-lmer(LarvaePerimeterum~ParentTrt*JarTrt+ (1|FemaleID)+(1|MaleID)+(1|CrossID)+(1|JarSeatable)+(1|JarID)+(1|TrtTankID)+(1|AccTankID), data=LarvaeDat)## boundary (singular) fit: see ?isSingularsteppedlarvperimFE<- step(larvperimFE)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00256065 (tol = 0.002, component 1)## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00531912 (tol = 0.002, component 1)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00469702 (tol = 0.002, component 1)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00531912 (tol = 0.002, component 1)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00281445 (tol = 0.002, component 1)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.004105 (tol = 0.002, component 1)print(steppedlarvperimFE)## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 12 -9768.4 19561
## (1 | TrtTankID) 1 11 -9768.4 19559 0.000 1 0.9964895
## (1 | AccTankID) 2 10 -9768.6 19557 0.436 1 0.5092703
## (1 | CrossID) 3 9 -9769.9 19558 2.517 1 0.1125893
## (1 | FemaleID) 0 8 -9788.8 19594 37.821 1 7.753e-10 ***
## (1 | MaleID) 0 8 -9777.2 19570 14.579 1 0.0001344 ***
## (1 | JarSeatable) 0 8 -9773.2 19562 6.639 1 0.0099792 **
## (1 | JarID) 0 8 -9882.5 19781 225.181 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## ParentTrt:JarTrt 0 425.2 425.2 1 240.57 5.9465 0.01547 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Model found:
## LarvaePerimeterum ~ ParentTrt + JarTrt + (1 | FemaleID) + (1 |
## MaleID) + (1 | JarSeatable) + (1 | JarID) + ParentTrt:JarTrt#final model chosen by the step function: LarvaePerimeterum~ParentTrt*JarTrt+ (1|FemaleID)+(1|MaleID)+(1|JarSeatable)+(1|JarID)
larvperimfinFE<- lmer(LarvaePerimeterum~ParentTrt*JarTrt+ (1|FemaleID)+(1|MaleID)+(1|JarSeatable)+(1|JarID), data=LarvaeDat)
#check the model to see that it meets assumptions
plot(larvperimfinFE) #seems to meet assumption of linearity and homoscedasticity
qqnorm(resid(larvperimfinFE))
qqline(resid(larvperimfinFE)) #has relatively long tails, doesn't meet assumption of normality. 
acf(resid(larvperimfinFE)) #this looks fine to me.
summary(larvperimfinFE)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: LarvaePerimeterum ~ ParentTrt * JarTrt + (1 | FemaleID) + (1 |
## MaleID) + (1 | JarSeatable) + (1 | JarID)
## Data: LarvaeDat
##
## REML criterion at convergence: 19539.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.3497 -0.5441 0.0477 0.6244 5.2637
##
## Random effects:
## Groups Name Variance Std.Dev.
## JarID (Intercept) 16.991 4.122
## FemaleID (Intercept) 8.259 2.874
## MaleID (Intercept) 3.925 1.981
## JarSeatable (Intercept) 1.430 1.196
## Residual 71.503 8.456
## Number of obs: 2698, groups:
## JarID, 270; FemaleID, 15; MaleID, 11; JarSeatable, 3
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 236.3700 1.6243 16.9133 145.519 <2e-16 ***
## ParentTrt2800 0.6594 2.0937 19.2644 0.315 0.7562
## JarTrtExposed -36.3906 0.8198 240.5488 -44.392 <2e-16 ***
## ParentTrt2800:JarTrtExposed 2.9251 1.1995 240.5746 2.439 0.0155 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) PrT2800 JrTrtE
## PrntTrt2800 -0.636
## JarTrtExpsd -0.256 0.199
## PrT2800:JTE 0.175 -0.288 -0.683Visualize the perimeter data
ggplot(LarvaeDat,
aes(fill = WCa_pHDIC, y = LarvaePerimeterum, x = JarOmegaCalcite)) +
geom_point(aes(shape= Block, colour=WCa_pHDIC)) +
ggtitle("Larvae area by larvae treatment")+
ylab("Larvae Perimeter (um)")+
theme_classic()
meanperim<- ddply(LarvaeDat, .(JarTrt, ParentTrt), numcolwise(mean, na.rm=T))
seperim<- ddply(LarvaeDat, .(JarTrt, ParentTrt), numcolwise(se, na.rm=T))
ggplot(data = meanperim, aes(x = ParentTrt, y = LarvaePerimeterum, fill=JarTrt)) +
geom_bar(stat="identity", aes(fill=JarTrt), position="dodge")+
labs(x = "Parent Treatment", y = "Larvae Perimeter (um)") +
scale_fill_manual(name = "Larvae Treatment", labels = c("Control", "OA Exposed"), values = c("skyblue2", "red2")) +
geom_errorbar(aes(ymin= LarvaePerimeterum-seperim$LarvaePerimeterum, ymax=LarvaePerimeterum+seperim$LarvaePerimeterum),width=0.2, position=position_dodge(.9))+
theme_classic()
Analyze survival data
surv1<- lmer(RatSurv~WCa_pHDIC+ (1|FemaleID)+(1|MaleID)+(1|AccTankID), data=SurvWideDat)## boundary (singular) fit: see ?isSingularsurvstep<- step(surv1)## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingularprint(survstep)## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 6 24.955 -37.911
## (1 | MaleID) 1 5 24.955 -39.911 0.0000 1 1.0000000
## (1 | AccTankID) 2 4 24.312 -40.624 1.2868 1 0.2566362
## (1 | FemaleID) 0 3 18.506 -31.012 11.6120 1 0.0006553 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## WCa_pHDIC 1 0.00030233 0.00030233 1 13 0.0289 0.8675
##
## Model found:
## RatSurv ~ (1 | FemaleID)#final model chosen by step function: RatSurv~ (1|FemaleID)
survfin<- lmer(RatSurv~ 1+ (1|FemaleID), data=SurvWideDat)
plot(survfin) #this looks fine to me
qqnorm(resid(survfin))
qqline(resid(survfin)) #I think this is fine
acf(resid(survfin))
summary(survfin)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: RatSurv ~ 1 + (1 | FemaleID)
## Data: SurvWideDat
##
## REML criterion at convergence: -51.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.37540 -0.64186 0.04267 0.52410 1.83497
##
## Random effects:
## Groups Name Variance Std.Dev.
## FemaleID (Intercept) 0.01122 0.1059
## Residual 0.01044 0.1022
## Number of obs: 45, groups: FemaleID, 15
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.86470 0.03131 14.00000 27.62 1.3e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#female ID and not parent treatment significantly affects survival
#I don't think a glmer is necessary this seems to meet the assumptions of lmer
survprop<- lmer(SurvChange~WCa_pHDIC+ (1|FemaleID)+(1|MaleID)+(1|AccTankID), data=SurvWideDat)## boundary (singular) fit: see ?isSingularsurvstepprop<- step(survprop)## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingularprint(survstepprop)## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 6 24.955 -37.911
## (1 | MaleID) 1 5 24.955 -39.911 0.0000 1 1.0000000
## (1 | AccTankID) 2 4 24.312 -40.624 1.2868 1 0.2566362
## (1 | FemaleID) 0 3 18.506 -31.012 11.6120 1 0.0006553 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## WCa_pHDIC 1 0.00030233 0.00030233 1 13 0.0289 0.8675
##
## Model found:
## SurvChange ~ (1 | FemaleID)#final model chosen by step function: SurvChange~ (1|FemaleID)
survfinprop<- lmer(SurvChange~ 1+ (1|FemaleID), data=SurvWideDat)
plot(survfinprop) #this looks fine to me
qqnorm(resid(survfinprop))
qqline(resid(survfinprop)) #I think this is fine
acf(resid(survfinprop))
summary(survfinprop)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: SurvChange ~ 1 + (1 | FemaleID)
## Data: SurvWideDat
##
## REML criterion at convergence: -51.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.83497 -0.52410 -0.04267 0.64186 2.37540
##
## Random effects:
## Groups Name Variance Std.Dev.
## FemaleID (Intercept) 0.01122 0.1059
## Residual 0.01044 0.1022
## Number of obs: 45, groups: FemaleID, 15
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.13530 0.03131 14.00000 4.322 0.000703 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#same result as survival proportion. Female ID is significant random effect but not parent treatmentAnalyze survival data with parent as a fixed effect
surv1FE<- lmer(RatSurv~ParentTrt+ (1|FemaleID)+(1|MaleID)+(1|AccTankID)+(1|TrtTankID), data=SurvWideDat)## boundary (singular) fit: see ?isSingularsurvstepFE<- step(surv1FE)## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingularprint(survstepFE)## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 7 24.769 -35.539
## (1 | MaleID) 1 6 24.769 -37.539 0.0000 1 1.0000000
## (1 | TrtTankID) 2 5 24.769 -39.539 0.0000 1 0.9999960
## (1 | AccTankID) 3 4 24.114 -40.228 1.3111 1 0.2522035
## (1 | FemaleID) 0 3 18.300 -30.601 11.6267 1 0.0006501 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## ParentTrt 1 0.0001736 0.0001736 1 13 0.0166 0.8994
##
## Model found:
## RatSurv ~ (1 | FemaleID)#final model chosen by step function: RatSurv~ (1|FemaleID), which is the same as survfin from above
survpropFE<- lmer(SurvChange~ParentTrt+ (1|FemaleID)+(1|MaleID)+(1|AccTankID)+(1|TrtTankID), data=SurvWideDat)## boundary (singular) fit: see ?isSingularsurvsteppropFE<- step(survpropFE)## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingularprint(survsteppropFE)## Backward reduced random-effect table:
##
## Eliminated npar logLik AIC LRT Df Pr(>Chisq)
## <none> 7 24.769 -35.539
## (1 | MaleID) 1 6 24.769 -37.539 0.0000 1 1.0000000
## (1 | TrtTankID) 2 5 24.769 -39.539 0.0000 1 0.9999960
## (1 | AccTankID) 3 4 24.114 -40.228 1.3111 1 0.2522035
## (1 | FemaleID) 0 3 18.300 -30.601 11.6267 1 0.0006501 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite
##
## Eliminated Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## ParentTrt 1 0.0001736 0.0001736 1 13 0.0166 0.8994
##
## Model found:
## SurvChange ~ (1 | FemaleID)#final model chosen by step function: SurvChange~ (1|FemaleID), which is the same as survfinpropVisualize larvae survival
par(mfcol=c(1,1))
boxplot(RatSurv~as.factor(ParentTrt), data=SurvWideDat)
boxplot(SurvChange~as.factor(ParentTrt), data= SurvWideDat)
#make reaction norms for survival color will be parent treatment. X will be jar treatment. Each cross will be its own line
ggplot(SurvDat,aes(x=JarTrt, y=LarvaeSurvived, color = ParentTrt, group= CrossID)) +
geom_point(aes(color=ParentTrt))+
geom_line(aes(color=ParentTrt, group=CrossID))+
geom_errorbar(aes(ymin= LarvaeSurvived-SELarvaeSurvived, ymax=LarvaeSurvived+SELarvaeSurvived),width=0.05, position=position_dodge(0))+
labs(x = "Larvae Treatment", y = "Number of Larvae Survived") +
scale_color_manual(values = c("skyblue2", "red2")) +
theme_classic()
ggplot(SurvWideDat,
aes(y = SurvChange, x = ParentTrt)) +
geom_point() +
ggtitle("Larvae survival by parent treatment")+
ylab("Proportion change in the number of larvae that survived")+
theme_classic()
meansurv<- ddply(SurvWideDat, .(ParentTrt), numcolwise(mean, na.rm=T))
sesurv<- ddply(SurvWideDat, .(ParentTrt), numcolwise(se, na.rm=T))
#make a shorter name for the CrossIDs in SurvWideDat
SurvWideDat$Cross<- substr(SurvWideDat$CrossID, 1, 9)
ggplot(data = SurvWideDat, aes(x = Cross, y = RatSurv, fill=ParentTrt)) +
geom_bar(stat="identity", position="dodge")+
labs(x = "Cross ID", y = "Proportion of Larvae that Survived") +
scale_fill_manual(name = "Parent Treatment", labels = c("Control", "OA Exposed"), values = c("skyblue2", "red2")) +
theme(axis.text.x=element_text(angle=90))
ggplot(data = meansurv, aes(x = ParentTrt, y = RatSurv)) +
geom_bar(stat="identity", position="dodge")+
labs(x = "Parent Treatment", y = "Proportion of Larvae that Survived")+
geom_errorbar(aes(ymin= RatSurv-sesurv$RatSurv, ymax=RatSurv+sesurv$RatSurv),width=0.2, position=position_dodge(.9))+
theme_classic()
Cilia extrusion statistics
#Plot cilia extrusion by size to see if there is a relationship (see Waldbusser et al 2015)
#to do this Elise looked at a subset of Larvae for each jar treatment.
par(mfcol=c(1,1))
testcil<- read.csv("../data/LarvMorphTest.csv")
testcil$CilExtruded<- as.factor(testcil$CilExtruded)
testcilsub<- subset(testcil, CilExtruded !="")
testcilsub$Flag<- as.factor(testcilsub$Flag)
testcilsub<- subset(testcilsub, Flag !="TRUE")
#take a look at the data
ggplot(testcilsub,
aes(fill = JarTrt, y = MaxFeretDiamum, x = CilExtruded)) +
geom_point(aes( shape=JarTrt, colour=JarTrt)) +
ggtitle("Larvae cilia extrusion by size")+
ylab("Larvae Diameter (um)")+
facet_grid(~JarTrt)+
theme_classic()
#test for a relationship between size and cilia extruded within each jar treatment type.
exp<- subset(testcilsub, JarTrt=="Exposed")
con<- subset(testcilsub, JarTrt=="Control")
expcil1<-glm(CilExtruded~MaxFeretDiamum, data=exp, family=binomial)
summary(expcil1) #no significant relationship between size and cilia extrusion##
## Call:
## glm(formula = CilExtruded ~ MaxFeretDiamum, family = binomial,
## data = exp)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.7172 0.7201 0.7444 0.7527 0.7725
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.566735 3.046585 0.186 0.852
## MaxFeretDiamum 0.008725 0.047015 0.186 0.853
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 182.22 on 163 degrees of freedom
## Residual deviance: 182.18 on 162 degrees of freedom
## AIC: 186.18
##
## Number of Fisher Scoring iterations: 4boxplot(MaxFeretDiamum~CilExtruded, data=exp)
#now test the control jars for relationship between size and cilia extrusion.
concil1<-glm(CilExtruded~MaxFeretDiamum, data=con, family=binomial)
summary(concil1) #looks like there is a significant relationship between size and cilia extrusion##
## Call:
## glm(formula = CilExtruded ~ MaxFeretDiamum, family = binomial,
## data = con)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.9324 -0.3705 -0.3098 -0.2538 2.8520
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 14.8599 7.5190 1.976 0.0481 *
## MaxFeretDiamum -0.2376 0.1026 -2.316 0.0206 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 82.805 on 179 degrees of freedom
## Residual deviance: 77.639 on 178 degrees of freedom
## AIC: 81.639
##
## Number of Fisher Scoring iterations: 6boxplot(MaxFeretDiamum~CilExtruded, data=con)
Export Figures
plots.dir.path<- list.files(tempdir(), pattern="rs-graphics", full.names= TRUE)
plots.png.paths<- list.files(plots.dir.path, pattern=".png", full.names=TRUE)
file.copy(from=plots.png.paths, to="../results")## logical(0)Look at weight per larvae
#we can check these data in a few ways: 1) look to see if replicate jars have similar weights; 2) check to see if first and second filters measured come up with similar values for weight per larva
#Start with 1)
#start by using all data:
boxplot(MeanWtPerLarvae~CrossID*JarTrt, data=LarvByJarDat)
#looks like some of the crosses have high variance. could we just eliminate those?
meanwtall<- ddply(LarvByJarDat, .(CrossID,ParentTrt, JarTrt, FemaleID, MaleID), numcolwise(mean, na.rm=T))
sewtall<- ddply(LarvByJarDat, .(CrossID, ParentTrt, JarTrt, FemaleID, MaleID), numcolwise(se, na.rm=T))
ggplot(data = meanwtall, aes(x = CrossID, y = MeanWtPerLarvae)) +
geom_errorbar(aes(ymin= MeanWtPerLarvae-sewtall$MeanWtPerLarvae, ymax=MeanWtPerLarvae+sewtall$MeanWtPerLarvae),width=0.2)+
geom_point(aes(colour=JarTrt)) +
labs(x = "Cross ID", y = "Mean Weight Per Larvae ug") +
theme_classic()
#let's see if removing filters with fewer than 100 individuals helps this out and those with crystals.
subdat<- LarvByJarDat
subdat$MeanWtPerLarvae<- ifelse(subdat$F1Crystals=="FALSE", paste(subdat$MeanWtPerLarvae), paste("NA"))
subdat$MeanWtPerLarvae<- ifelse(subdat$F1LarvaeCount>100, paste(subdat$MeanWtPerLarvae), paste("NA"))
subdat$MeanWtPerLarvae<- ifelse(subdat$JarID!="J004_B1", paste(subdat$MeanWtPerLarvae), paste("NA"))#this is the only F2 filter after subsetting that has <100 larvae
subdat$JarSeatable<- as.factor(subdat$JarSeatable)
subdat$MeanWtPerLarvae<- as.numeric(subdat$MeanWtPerLarvae)## Warning: NAs introduced by coercionsubdat$JarTrt<- as.factor(subdat$JarTrt)
subdat$LarvDensity<- subdat$MeanWtPerLarvae/subdat$LarvaeAreaum2
boxplot(MeanWtPerLarvae~CrossID*JarTrt, data=subdat)
#looks like some of the crosses have high variance. could we just eliminate those?
submeanwtall<- ddply(subdat, .(Block, ParentTrt, JarTrt, CrossID, FemaleID, MaleID), numcolwise(mean, na.rm=T))
subsewtall<- ddply(subdat, .(Block, ParentTrt, JarTrt, CrossID, FemaleID, MaleID), numcolwise(se, na.rm=T))
ggplot(data = submeanwtall, aes(x = CrossID, y = MeanWtPerLarvae)) +
geom_errorbar(aes(ymin= MeanWtPerLarvae-subsewtall$MeanWtPerLarvae, ymax=MeanWtPerLarvae+subsewtall$MeanWtPerLarvae),width=0.2)+
geom_point(aes(colour=JarTrt)) +
labs(x = "Cross ID", y = "Mean Weight Per Larvae ug") +
theme(axis.text.x = element_text(angle = 90))## Warning: Removed 4 rows containing missing values (geom_point).
#there are some crosses that only have a single point. Those should be removed. Do this by looking at se
#first join SE to the dataset
colnames(subsewtall)<- paste("SE", colnames(subsewtall), sep="")
#merge the standard error data to the SurvDat dataframe
SubWt<- merge(submeanwtall, subsewtall, by.x= c("CrossID", "JarTrt"), by.y= c("SECrossID", "SEJarTrt"))
relcrosses<- SubWt
relcrosses$MeanWtPerLarvae<- ifelse(relcrosses$SEMeanWtPerLarvae>0, paste(relcrosses$MeanWtPerLarvae), paste("NA"))
relcrosses$MeanWtPerLarvae<- as.numeric(relcrosses$MeanWtPerLarvae)
ggplot(data = relcrosses, aes(x = CrossID, y = MeanWtPerLarvae)) +
geom_errorbar(aes(ymin= MeanWtPerLarvae-SEMeanWtPerLarvae, ymax=MeanWtPerLarvae+SEMeanWtPerLarvae),width=0.2)+
geom_point(aes(colour=JarTrt)) +
labs(x = "Cross ID", y = "Mean Weight Per Larvae ug") +
theme(axis.text.x = element_text(angle = 90))## Warning: Removed 22 rows containing missing values (geom_point).
boxplot(SEMeanWtPerLarvae~JarTrt, data= relcrosses)
#I don't know how to choose which ones to include. For now, let's say if the SE is greater than mean(MeanWtPerLarvae) then I will remove it
MeanWtVal<- mean(relcrosses$MeanWtPerLarvae, na.rm=TRUE)
relcrosses$MeanWtPerLarvae<- ifelse(relcrosses$SEMeanWtPerLarvae<= MeanWtVal, paste(relcrosses$MeanWtPerLarvae), paste("NA"))
relcrosses$MeanWtPerLarvae<- as.numeric(relcrosses$MeanWtPerLarvae)## Warning: NAs introduced by coercionboxplot(MeanWtPerLarvae~ParentTrt*as.factor(JarTrt), data=relcrosses)
ggplot(data = relcrosses, aes(x = CrossID, y = LarvDensity)) +
geom_errorbar(aes(ymin= LarvDensity-SELarvDensity, ymax=LarvDensity+SELarvDensity),width=0.2)+
geom_point(aes(colour=JarTrt)) +
labs(x = "Cross ID", y = "Mean Density Per Larvae ug/um2") +
theme(axis.text.x = element_text(angle = 90))## Warning: Removed 4 rows containing missing values (geom_point).
#now try to create a linear model
Density1<- lm(LarvDensity~ParentTrt*as.factor(JarTrt), data=relcrosses)
summary(Density1)##
## Call:
## lm(formula = LarvDensity ~ ParentTrt * as.factor(JarTrt), data = relcrosses)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.607e-10 -2.152e-10 -5.552e-11 1.227e-10 1.144e-09
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.929e-10 7.358e-11 6.698 2.45e-09
## ParentTrt2800 -1.336e-11 1.091e-10 -0.122 0.9029
## as.factor(JarTrt)Exposed 2.195e-10 1.064e-10 2.063 0.0423
## ParentTrt2800:as.factor(JarTrt)Exposed 1.515e-10 1.559e-10 0.972 0.3340
##
## (Intercept) ***
## ParentTrt2800
## as.factor(JarTrt)Exposed *
## ParentTrt2800:as.factor(JarTrt)Exposed
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.605e-10 on 82 degrees of freedom
## (4 observations deleted due to missingness)
## Multiple R-squared: 0.1598, Adjusted R-squared: 0.129
## F-statistic: 5.198 on 3 and 82 DF, p-value: 0.002452#but from the plot perhaps something is relevant for the difference going from control to exposed conditions.
#use spread function
relcrossestowide<- unite(relcrosses, Value, c(MeanWtPerLarvae, SEMeanWtPerLarvae, LarvDensity, SELarvDensity), sep="_")
relcrossestowide<- subset(relcrossestowide, select=c("CrossID", "Value", "JarTrt", "ParentTrt", "MaleID", "FemaleID"))
#use spread function in tidyr to make the column in long form vs. wide form
relwide<- spread(relcrossestowide, JarTrt, value= Value)
#now split the mean and SE apart again
relwide<- separate(relwide, Control, into= c("MeanWtPerLarvaeCon", "SEMeanWtPerLarvaeCon", "LarvDensityCon", "SELarvDensityCon"), sep="_")
relwide<- separate(relwide, Exposed, into= c("MeanWtPerLarvaeExp", "SEMeanWtPerLarvaeExp", "LarvDensityExp", "SELarvDensityExp"), sep="_")
#make the means and ses numeric
relwide$MeanWtPerLarvaeCon<- as.numeric(relwide$MeanWtPerLarvaeCon)## Warning: NAs introduced by coercionrelwide$MeanWtPerLarvaeExp<- as.numeric(relwide$MeanWtPerLarvaeExp)## Warning: NAs introduced by coercionrelwide$SEMeanWtPerLarvaeCon<- as.numeric(relwide$SEMeanWtPerLarvaeCon)## Warning: NAs introduced by coercionrelwide$SEMeanWtPerLarvaeExp<- as.numeric(relwide$SEMeanWtPerLarvaeExp)## Warning: NAs introduced by coercionrelwide$LarvDensityCon<- as.numeric(relwide$LarvDensityCon)
relwide$SELarvDensityCon<- as.numeric(relwide$SELarvDensityCon)## Warning: NAs introduced by coercionrelwide$LarvDensityExp<- as.numeric(relwide$LarvDensityExp)
relwide$SELarvDensityExp<- as.numeric(relwide$SELarvDensityExp)## Warning: NAs introduced by coercion#get the difference from control to OA
relwide$WtDiff<- (relwide$MeanWtPerLarvaeExp-relwide$MeanWtPerLarvaeCon)
plot(WtDiff~ParentTrt, data=relwide)
boxplot(WtDiff~FemaleID, data=relwide)
abline(a=0, b=0)
relwide$DensDiff<- (relwide$LarvDensityExp-relwide$LarvDensityCon)
plot(DensDiff~ParentTrt, data=relwide)
#now try a linear model.
Dens1<- lm(DensDiff~ParentTrt, data=relwide)
summary(Dens1)##
## Call:
## lm(formula = DensDiff ~ ParentTrt, data = relwide)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.901e-10 -1.801e-10 -1.457e-11 1.668e-10 7.025e-10
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.007e-10 6.006e-11 3.342 0.00185 **
## ParentTrt2800 1.184e-10 8.822e-11 1.342 0.18751
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.817e-10 on 39 degrees of freedom
## (4 observations deleted due to missingness)
## Multiple R-squared: 0.04411, Adjusted R-squared: 0.0196
## F-statistic: 1.8 on 1 and 39 DF, p-value: 0.1875#Look at the difference between F1 and F2 measurements, so subset the data for only those with F2 measurements
#get the mean of weight per larva
mean(LarvByJarDat$MeanWtPerLarvae)## [1] NAF1F2<- subset(LarvByJarDat, F2WtPerLarvae >0)
F1F2<- subset(F1F2, F1WtPerLarvae>0)
F1F2$JarTrt<- as.factor(F1F2$JarTrt)
F1F2$F1FilterID<- as.factor(F1F2$F1FilterID)
F1F2$F2FilterID<- as.factor(F1F2$F2FilterID)
F1F2$ParentTrt<- as.factor(F1F2$ParentTrt)
F1F2$CrossID<- as.factor(F1F2$CrossID)
F1F2$TrtTankID<- as.factor(F1F2$TrtTankID)
F1F2$JarSeatable<- as.factor(F1F2$JarSeatable)
F1F2$MaleID<- as.factor(F1F2$MaleID)
F1F2$FemaleID<- as.factor(F1F2$FemaleID)
F1F2$JarID<-droplevels(F1F2)$JarID
F1F2<- subset(F1F2, select=c("JarID", "F1WtPerLarvae", "F1LarvaeCount", "F2WtPerLarvae", "F2LarvaeCount", "F1FilterID", "F2FilterID", "JarTrt", "ParentTrt", "CrossID", "TrtTankID", "JarSeatable", "MaleID", "FemaleID", "F1Crystals", "MeanWtPerLarvae"))
F1F2$LarvWtDiff<- abs(F1F2$F1WtPerLarvae-F1F2$F2WtPerLarvae)
F1F2$LarvWtDiffProp<- (F1F2$LarvWtDiff)/ (F1F2$F2WtPerLarvae)
F1F2$TotalLarvCount<- F1F2$F1LarvaeCount+F1F2$F2LarvaeCount
F1F2$LarvaeCountDiff<- abs(F1F2$F1LarvaeCount- F1F2$F2WtPerLarvae)
F1F2sub<- subset(F1F2, F1LarvaeCount>100)
F1F2sub<- subset(F1F2sub, F2LarvaeCount>100)
F1F2sub<- subset(F1F2sub, F1Crystals=="FALSE")
F1F2sub$JarID<-droplevels(F1F2sub)$JarID
#Now look at the differences.In weights
mean(F1F2sub$F1WtPerLarvae)## [1] 2.213654e-06mean(F1F2sub$F2WtPerLarvae)## [1] 1.594731e-06plot(LarvWtDiff~TotalLarvCount, data=F1F2sub, col="black", pch=19)
abline(a= 1.59e-06, b=0, col="blue")
abline(a=(2* 1.59e-06), b=0, col="red")
plot(F1LarvaeCount~LarvWtDiffProp, data=F1F2sub, pch=19)
points(F2LarvaeCount~LarvWtDiffProp, data=F1F2sub, pch=19, col="orange")
#maybe I should just remove those with larvae counts <200
F200sub<- subset(F1F2, F1LarvaeCount>200)
F200sub<- subset(F200sub, F2LarvaeCount>200)
F200sub$JarID<-droplevels(F200sub)$JarID
mean(F200sub$F1WtPerLarvae)## [1] 1.91e-06mean(F200sub$F2WtPerLarvae)## [1] 1.256176e-06plot(LarvWtDiff~TotalLarvCount, data=F200sub, col="black", pch=19) 
plot(F1LarvaeCount~LarvWtDiffProp, data=F200sub, pch=19)
points(F2LarvaeCount~LarvWtDiffProp, data=F200sub, pch=19, col="orange")
plot(F1WtPerLarvae~JarID, data=F200sub)
points(F2WtPerLarvae~JarID, data=F200sub)
points(LarvWtDiff~JarID, data=F200sub, col="red", pch=19) 
boxplot(LarvWtDiffProp~ParentTrt*JarTrt, data=F1F2sub)
#I'm going to just look at the jars that had proportion difference <1
Less1<- subset(F1F2, LarvWtDiffProp<1)
boxplot(MeanWtPerLarvae~ParentTrt*JarTrt, data=Less1)
boxplot(LarvWtDiffProp~ParentTrt*JarTrt, data=Less1)
Wtdiff<- lm(LarvWtDiff~TotalLarvCount, data=F1F2)
par(mfcol=c(2,2))
plot(Wtdiff)
par(mfcol=c(1,1))
summary(Wtdiff)##
## Call:
## lm(formula = LarvWtDiff ~ TotalLarvCount, data = F1F2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.642e-06 -4.206e-06 -2.710e-06 1.996e-06 3.599e-05
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.669e-06 2.454e-06 3.94 0.000347 ***
## TotalLarvCount -6.907e-09 3.153e-09 -2.19 0.034877 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.982e-06 on 37 degrees of freedom
## Multiple R-squared: 0.1148, Adjusted R-squared: 0.09086
## F-statistic: 4.798 on 1 and 37 DF, p-value: 0.03488#Use the Wt diff model to decide which filters to include
#first get the average larva weight overall:
mean(F1F2$MeanWtPerLarvae)## [1] 3.874513e-06boxplot(F1WtPerLarvae~ParentTrt+JarTrt, data=Calc)## Error in eval(m$data, parent.frame()): object 'Calc' not found#I don't know whether or not to trust these data. Let's look at F1 and F2 together
Calcmean<- ddply(Calc, .(ParentTrt, JarTrt), numcolwise(mean, na.rm=TRUE))## Error in empty(.data): object 'Calc' not foundCalcse<- ddply(Calc, .(ParentTrt, JarTrt), numcolwise(se, na.rm=TRUE))## Error in empty(.data): object 'Calc' not foundggplot(data = Calcmean, aes(x = ParentTrt , y = F1WtPerLarvae, fill=JarTrt)) +
geom_bar(stat="identity", aes(fill=JarTrt), position="dodge")+
labs(x = "Parent Treatment", y = "Weight per Larvae") +
scale_fill_manual(name = "Larvae Treatment", labels = c("Control", "OA Exposed"), values = c("skyblue2", "red2")) +
geom_errorbar(aes(ymin= F1WtPerLarvae-Calcse$F1WtPerLarvae, ymax=F1WtPerLarvae+Calcse$F1WtPerLarvae),width=0.2, position=position_dodge(.9))+
theme_classic()## Error in ggplot(data = Calcmean, aes(x = ParentTrt, y = F1WtPerLarvae, : object 'Calcmean' not foundcalcexp<- subset(Calc, JarTrt=="Exposed")## Error in subset(Calc, JarTrt == "Exposed"): object 'Calc' not foundcalcexp$JarTrt<- as.factor(calcexp$JarTrt)## Error in is.factor(x): object 'calcexp' not foundcalccon<- subset(Calc, JarTrt=="Control")## Error in subset(Calc, JarTrt == "Control"): object 'Calc' not foundcalccon$JarTrt<- as.factor(calccon$JarTrt)## Error in is.factor(x): object 'calccon' not foundplot(F1WtPerLarvae~F1LarvaeCount, data=calccon, col="blue", pch=19)## Error in eval(m$data, eframe): object 'calccon' not foundpoints(F1WtPerLarvae~F1LarvaeCount, data=calcexp, col="red", pch=19)## Error in eval(m$data, eframe): object 'calcexp' not foundsummary(lm(F1WtPerLarvae~F1LarvaeCount+JarTrt+ParentTrt, data=Calc))## Error in is.data.frame(data): object 'Calc' not found